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ABSTRACT 

A  genera]  analysis  applicable  to  a  variety  of  free  mixing  flows  of  engineering  interest 
is  described.  An  empirical  relationship  between  the  turbulent  shear  stress  and  the  turbulent 
kinetic  energy  is  used.  Solution  of  the  turbulent  kinetic  energy  equation  as  one  of  the 
governing  equations  yields  the  turbulent  shear  stress  at  all  points  in  the  flow  field.  Algebraic 
relationships  are  used  to  obtain  the  length  scales  required  to  close  the  turbulent  kinetic 
energy  equation.  The  computer  program  developed  to  provide  the  numerical  framework 
for  this  analysis  is  described,  and  a  FORTRAN  listing  and  description  of  required  inputs 
is  included.  To  demonstrate  the  capabilities  of  the  analysis,  a  number  of  experimental 
free  mixing  flows  have  been  calculated,  and  the  results  of  these  calculations  are  discussed. 
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NOMENCLATURE 

a  Coefficient  in  general  transport  equation 

ai  Ratio  of  turbulent  shear  stress  to  turbulent  kinetic  energy 

a2  Turbulent  kinetic  energy  dissipation  coefficient 

b  Coefficient  in  general  transport  equation 

C  Jet  species  concentration 

c  Coefficient  in  general  transport  equation 

Ci  Density-ratio  correction  to  ai 

d  Coefficient  in  general  transport  equation 

H  Total  enthalpy 

k  Turbulent  kinetic  energy  (TKE) 

kp  Constant  in  Prandtl  eddy  viscosity  model 

8k  TKE  dissipation  length  scale 

ra"  Entrained  mass  flux/unit  area 

P  Pressure 

Pr-j  Turbulent  Prandtl  number 
Pr*  Turbulent  Prandtl  number  for  TKE 
R  Gas  constant 
r  Radial  coordinate 

rj  Radius  of  inner  edge  of  mixing  region 

Rt  Turbulent  Reynolds  number 
Sex  Turbulent  Schmidt  number 
Tt  Total  temperature 
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u  Axial  mean  velocity  component 

iT  Velocity  desired  at  a  flow  field  edge,  entrainment  calculation 
uc  Mean  velocity  on  jet  centerline  at  any  x 
ue  Mean  velocity  in  surrounding  stream 
Uj  Mean  velocity  on  jet  centerline  at  nozzle  exit 
v  Lateral  mean  velocity  component 
x  Axial  coordinate 
y  Lateral  coordinate 
Greek  Symbols 

a  Parameter;  =  0  for  planar  flow,  =  1  for  axisymmetric  flow 
Au  Characteristic  velocity  difference 
e  Turbulent  eddy  viscosity 

em  Turbulent  eddy  viscosity  at  a  characteristic  point 
8  Angle  between  flow  and  x-coordinate 

p  Fluid  density 

pci  Density  of  surrounding  stream  at  nozzle  exit 
Pji  Density  of  jet  at  nozzle  exit 
a  Two-dimensional  shear-layer  spread  parameter 
a o  Value  of  a  for  constant-density,  single-stream  shear  layer 
Transport  coefficient  ratio  for  general  variable  <p 
t  Turbulent  shear  stress 
<p  General  transportable  variable 
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SECTION  1 
INTRODUCTION 

Free  turbulent  mixing  phenomena  are  of  fundamental  importance  in  many  flow 
systems  of  engineering  interest.  Applications  in  which  free  turbulent  mixing  phenomena 
occur  include  aircraft  and  missile  plumes,  jet  pumps,  wakes,  and  combustors.  The  correct 
prediction  of  such  flows  is  a  basic  requirement  for  the  solution  of  many  more  complex 
problems,  such  as  those  involving  chemical  reactions  in  a  flowing  system  and  prediction 
of  exhaust  jet  noise.  The  general  theory  of  free  turbulent  mixing  described  in  this  report 
has  been  developed  to  provide  an  accurate  prediction  of  both  the  overall  mixing  rate 
in  a  free  turbulent  flow  and  the  details  of  the  mixing  process,  such  as  the  profiles  of 
the  turbulent  shear  stress,  velocity,  total  enthalpy,  and  jet  species.  The  theory  is  applicable 
to  a  broad  range  of  free  turbulent  flows,  including  constant  and  variable  density  jets  and 
wakes;  the  required  expenditure  of  computer  time  is  small,  approximately  one  minute 
for  a  typical  case  on  an  IBM  370/155  computer. 

1.1  CAPABILITIES  AND  LIMITATIONS 

The  basic  approach  for  the  analysis  of  the  free  turbulent  mixing  problems  described 
in  this  report  is  developed  from  the  turbulent  kinetic  energy  (TKE)  method  described 
in  Ref.  1.  The  method  incorporates  several  of  the  improvements  which  were  discussed 
in  Ref.  2,  as  well  as  certain  modifications  which  have  been  incorporated  since  the  latter 
description  was  written.  These  modifications  allow  the  use  of  a  single,  consistent  model 
for  all  of  the  calculations  described  in  Ref.  2.  The  technique  allows  the  computation 
of  planar  and  axisymmetric  one-  and  two-stream  mixing  problems.  Referring  to  Fig.  1 
as  a  schematic  of  the  types  of  flow  field  involved,  any  combination  of  velocities  ue  and 
Uj  may  be  considered,  except  ue  or  uj  =  0.  Single-stream  flows  are  approximated  by  small 


Fig.  1  Schematic  of  Free  Mixing  Flow  Field 
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values  of  the  velocity  ratio,  typically  Ue/Uj  ^  0.01.  The  densities  of  the  two  streams  may 
be  different,  and  computations  have  been  made  successfully  for  values  of  the  density  ratio 
pe/pj  as  great  as  14.  Such  cases,  in  which  the  center  stream  density  is  considerably  lower 
than  the  outer  stream  density,  are  the  most  difficult  to  handle  because  of  numerical 
difficulties  associated  with  the  coordinate  transformation;  no  limitations  have  been  observed 
for  Pe/P]  «  1.  Transport  of  energy  and  mass,  as  well  as  momentum,  can  be  handled, 
and  nonunity  turbulent  Prandtl  and  Schmidt  numbers  can  be  assumed.  As  presently 
formulated,  the  turbulent  Prandtl  and  Schmidt  numbers  are  assumed  to  be  constant 
everywhere,  but  there  is  no  fundamental  reason  why  axial  and  radial  variations  of  these 
transport  coefficients  could  not  be  prescribed. 

The  technique  is  specialized  to  free  turbulent  processes,  i.e.,  the  mixing  field  is 
assumed  to  exist  in  a  constant-pressure  environment  remote  from  walls.  A  longitudinal 
pressure  gradient  may  be  incorporated  with  some  modifications  to  the  program,  but  it 
must  be  specified  a  priori.  Although  compressible  flows  can  be  handled,  they  must  be 
shock  free;  in  general,  this  limits-  the  computation  of  supersonic  flows  to  those  from 
correctly  expanded  nozzles.  Only  mathematically  parabolic  flows  of  the  boundary-layer 
type  are  considered.  In  general  this  limitation  implies  that  the  nozzle  lip  shown  in  Fig. 
1  is  thin;  otherwise  a  region  of  recirculating  "(and,  therefore,  mathematically  elliptic)  flow 
will  exist  just  downstream  of  the  lip.  In  practice,  this  is  not  a  serious  limitation  since 
the  assumption  of  a  region  of  small  but  positive  velocity  in  the  lip  region  will  normally 
suffice  for  accurate  prediction  of  the  mixing  field.  Finally,  the  approach  is  limited  to 
the  mixing  of  two  streams.  Again,  there  is  jno  fundamental  reason  why  more  than  two 
streams  could  not  be  handled.  However,  such  a  multi-stream  mixing  process  would  require 
a  careful  analysis  of  the  problem  of  merging  shear  layers  with  different  characteristic  length 
scales;  in  turn,  this  would  require  additional,  fairly  complex  programming  logic. 

1.2  A  BRIEF  HISTORY  OF  APPROACHES 

The  fundamental  problem  in  the  computation  of  a  turbulent  shear  flow  is  the 
specification  of  thd  turbulent  shear  stress,  t.  Attempts  to  solve  this  problem  have  occupied 
a  great  number  of  research  workers  for  nearly  a  century.  It  was  only  recently  with  the 
availability  of  high-speed  computers,  that  reliable  computations  of  any  but  the  simplest 
free  turbulent  flows  have  been  possible. 

Turbulent  shear  stress  models  are  discussed  at  some  length  in  Ref.  1 ,  so  only  a  general 
outline  will  be  given  here.  Such  models  can  be  divided  into  two  classes:  those  in  which 
the  turbulent  shear  stress  is  taken  to  be  a  function  of  the  local  mean  velocity  and  a 
characteristic  mean-flow  length  scale,  and  those  in  which  a  transport  equation  is  solved 
to  obtain  the  turbulent  shear  stress  simultaneously  with  the  mean-flow  equations  describing 
the  problem. 

The  first  class  of  models  includes  the  mixing  length  (Ref.  and  eddv.  viscosity  (Ref. 
4)  models  of  Prandtl,  as  well  as  the  vorticity  transport  theory  of  Taylor  (Ref.  5) 
and  the  "inductive"  approach  of  Reichardt  (Ref.  6).  More  recently,  the  "displacement 
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thickness"  model  has  been  proposed  by  Schetz  (Ref.  7).  All  of  these  models  share  a 
common  fault.  Although  all  of  them  are  capable  of  computing  one  flow  or  a  few  flows 
well  (generally  those  flows  which  were  used  to  evaluate  the  inevitable  empirical  constants 
that  exist  in  the  models),  they  cannot  be  used  with  confidence  in  any  other  flow  (Ref. 
1).  This  failure  of  eddy  viscosity  models  has  tended  to  retard  the  use  of  theoretical  models 
for  prediction  of  the  free  turbulent  flows  encountered  in  engineering  practice. 

The  second  class  of  models  includes  the  eddy-viscosity  transport  equation  developed 
by  Nee  and  Kovasznay  (Ref.  8)  as  well  as  several  turbulent  kinetic  energy  models.  The 
models  are  based  on  expressions  for  the  turbulent  shear  stress  which  either  relate  it  to  the 
turbulent  kinetic  energy  (such  as  the  present  work)  or  involve  a  solution  of  a 
turbulent  shear  stress  transport  equation.  In  recent  years  a  number  of  TKE  models  have 
been  developed.  These  models  can  be  arranged  in  a  hierarchy,  ranked  by  order  of 
complexity,  which  ranges  from  the  solution  of  one  additional  transport  equation,  for  the 
turbulent  kinetic  energy,  to  the  solution  of  transport  equations  for  the  turbulent  shear 
stress,  all  components  of  the  turbulent  kinetic  energy,  and  a  characteristic  turbulent  length 
scale.  The  present  model  fits  into  the  first  level  of  this  hierarchy,  involving  the  solution 
of  the  TKE  transport  equation  along  with  an  algebraic  formulation  for  the  turbulent  length 
scale  which  enters  the  problem.  The  model  proposed  by  Rodi  and  Spalding  (Ref.  9)  is 
at  the  next  level  of  complexity,  involving  a  TKE  equation  and  a  length  scale  equation. 
At  the  most  complex  level,  the  approach  of  Donaldson  and  Rosenbaum  (Ref.  10)  involves 
the  solution  of  five  equations  in  addition  to  the  mean  momentum  equation  (common  to 
all  of  these  approaches)  for  even  the  simplest  problem. 

At  no  point  in  the  hierarchy  is  the  need  for  empirical  input  eliminated.  There  appears 
to  be  an  almost  religious  belief  professed  by  researchers  involved  in  turbulence  that  as 
the  number  of  equations  increases,  so  does  the  profundity  of  the  results;  but  the  fact 
remains  that  as  the  number  of  equations  increases,  the  number  of  empirical  parameters 
increases  also,  and  the  physical  basis  for  the  necessary  empirical  input  becomes  more 
obscure.  In  the  present  model  the  empiricism  is  used  to  make  the  TKE  equation  soluble, 
and  since  the  turbulent  length  scales  appear  to  be  directly  related  to  the  shear-layer  width 
for  most  free  turbulent  flows,  algebraic  formulae  are  adopted  to  relate  the  characteristic 
length  scale  to  the  shear-layer  width.  The  results  of  the  1972  NASA-Langley  Conference 
on  Free  Turbulent  Shear  Flows  (Ref.  11)  showed  that  despite  the  relative  simplicity  of 
the  assumptions  used  in  this  model,  the  present  model  and  a  closely  related  one  (Ref. 
12)  were  two  of  the  most  versatile  methods  presented. 

1.3  LAYOUT  OF  THE  REPORT 

In  this  report  the  details  of  the  finite-difference  turbulent  kinetic  energy  (TKE)  model 
developed  by  the  author  will  be  described.  Following  this  description,  a  summary  of  the 
numerical  framework  used  to  solve  the  equations  of  motion  (based  on  the  technique 
developed  by  Patankar  and  Spalding  (Ref.  13))  will  be  presented.  Results  of  the  method 
as  they  applied  to  the  test  cases  from  the  1972  NASA-Langley  Conference  will  then  be 
described,  followed  by  a  user's  manual  for  the  computer  program  and  a  complete  program 
listing. 
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SECTION  II 
ANALYSIS 

2.1  THE  GOVERNING  EQUATIONS 


The  basic  equations  to  be  solved  in  any  analysis  of  a  free  turbulent  flow  (such  as 
that  shown  schematically  in  Fig.  1)  are  the  turbulent  transport  equations  for  mass, 
momentum,  energy,  and  species.  To  this  set,  the  present  analysis  adds  the  turbulent  kinetic 
energy  equation,  as  well  as  several  algebraic  relationships  between  the  dependent  variables. 
Because  all  free  turbulent  mixing  problems  to  be  considered  involve  flow  fields  in  which 
gradients  with  y  are  very  much  greater  than  those  with  x,  the  boundary-layer 
approximations  may  be  applied.  Thus  the  governing  Navier-Stokes  equations  are  reduced 
to  the  following  set  of  parabolic  equations: 

Continuity: 


<?u 

'"‘is; + 


-° 


Momentum: 


Species  Transport: 


+  pv 


<?u  1  d 

dy  ~  yt  dy 


(y°r)  - 


dp 

dx 


<?c  do  l  d  (pi  ya  <9C\ 

Pa  d*  +  pV  dy  ~  ya  dy  \  Sc  dy) 


where 


e  =  (r/p)/(5  u/(9y) 


(1) 


(2) 


(3) 


(4) 


■v 


Mean  Energy  Transport: 


pu 


5H 

5x 


+ 


]_d_ 
ya  dy 


(5) 


Turbulent  Kinetic  Energy: 


(6) 


For  planar  flow,  a  =  0,  while  for  axisymmetric  flow,  a  =  1. 
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In  writing  the  species,  mean  energy,  and  kinetic  energy  transport  equations,  gradient 
diffusion  has  been  assumed.  A  similar  form  for  the  diffusion  term  in  the  momentum 
equation  can  be  obtained  through  the  use  of  the  Boussinesq  hypothesis  (Eq.  (4));  thus. 


and  all  of  the  equations  may  be  written  in  the  standard  form 


It  must  be  emphasized  that  Eq.  (4)  serves  in  the  present  analysis  merely  to  provide 
a  relationship  with  whicli  to  evaluate  a  momentum  diffusion  coefficient.  Unlike  other 
analyses  of  turbulent  flows,  this  coefficient  is  not  directly  modeled  in  order  to  obtain 
the  turbulent  shear  stress,  but  it  serves  simply  as  a  formalism  which  allows  all  of  the 
equations  of  motion  to  be  written  in  the  form  of  Eq.  (8).  A  model  for  the  turbulent 
shear  stress  must  still  be  specified,  and  this  specification  is  the  subject  of  the  following 
section. 

2.2  SHEAR  STRESS  MODEL 

Equations  (1)  through  (6)  do  not  constitute  a  closed  set  of  equations  since  there 
is  no  relationship  yet  specified  for  the  turbulent  shear  stress,  r.-  In  this  work,  a  relationship 
proposed,  after  Nevzgljadov,  by  Dryden  (Ref.  14)  and  used  in  computations  of  the 
turbulent  boundary  layer  by  Bradshaw,  et  al.  (Ref.  15),  is  used: 

r  =  a]Pk  (9) 

Although  there  is  evidence  that  in  strong  shear  flows  the  parameter  ai  is  constant  with 
a  value  of  0.3  over  most  of  the  flow  (Ref.  16).  in  certain  regions  of  the  flow  field  it 
cannot  be  constant.  For  example,  in  the  region  of  the  centerline  in  an  axisynunetric  flow, 
the  turbulent  shear  stress  must  apprqach  zero,  yet  the  turbulent  kinetic  energy  does 
not.  Similarly,  in  the  two-dimensional  free  shear  layer,  large  variations  of  ai  are 
experimentally  observed  near  the  edges  of  the  flow.  For  these  reasons,  a 
semiempirical  set  of  profiles  describing  the  variation  of  aj  from  its  nominal  value  of  0.3 
has  been  adopted  for  this  work.  The  relationships  used  are  defined  in  Table  1  (Eqs.  (10) 
through  (13)). 

The  shape  of  the  a  profile  predicted  by  Eqs.  (12)  and  (13)  (Table  1)  agrees  quite 
well  with  experimental  data  for  the  variation  of  ai  in  incompressible  jets,  as  Fig.  2  shows. 
Equation  (10)  has  been  verified  by  the  accurate  prediction  of  experimental  velocity  profiles, 
and  Eq.  (11)  is  inserted  at  the  beginning  of  the  second  regime  as  a  simple  approximation 
to  the  actual  transition  between  Eq.  (10)  and  Eqs.  (12)  and  (13).  Little  experimental 
evidence  is  available  for  the  actual  shape  of  the  parameter  a]  through  this  transition  region. 


5 


AED  C-T  R-73-1  77 


TABLE  I 

PROFILE  RELATIONS  FOR  THE  PARAMETER  3l 


Flow  Regime 


Range 


Expression  for  a. 


Eq. 


Remarks 


0.05  <- 


u  -  u 

J  e 


<0.95 


2-D  Shear -Layer  First 
Regime  of  Jets  (Fig.  1) 


u  -  u 

0  < - —  <  0.05 

—  u  -  u 

J  e 
u  -  u 

0.  95  < - —  <1.0 

u  -  u 
J  e 


None 


(10) 


None 


Ratio  of  velocity  gradients  provides 
proper  sign 


e  held  constant  at  value  obtained  from 
Eqs.  (4)  and  (9>  at  (u  -  ue)/(uj  -  uc)  -  0.05, 
r  obtained  from  Eq.  (4) 

c  held  constant  at  value  obtained  from 
Eqs.  (4)  and  <9)  at  <u  -  ue)/(uj  -  u^)  i  0.95, 
t  obtained  from  Eq.  (4) 


Second  Regime  of 
A  xi  symmetric  and 
Planar  Jets  (Fig.  1), 
Wakes 


u  >  0.  9  u 
c  ] 


u  <0,9ut 
c  J 


uc  0.  9  Uj 
y  >  ^max 


<u) 


(12) 


(13) 


Ratio  of  velocity  gradients  provides 
proper  sign 

ym  is  the  value  of  lateral  coordinate 
at  which  [fiu/Ayl  is  a  maximum 


Ratio  of  velocity  gradients  provides 
proper  sign 


1.1 


Fig.  2  Comparison  of  Calculated  and  Experimental  Variation 
of  Ratio  of  Shear  Stress  to  Kinetic  Energy 
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2.3  MODELING  OF  THE  TURBULENT  KINETIC  ENERGY  EQUATION 
2.3.1  The  Convection  and  Production  Terms 


The  turbulent  kinetic  energy  (TKE)  equation  as  written  in  Section  2.1  (Eq.  (6)) 
incorporates  several  assumptions.  As  with  .any  equation  for  a  transportable  scalar,  this 
equation  can  be  written 

convection  =  diffusion  +  production  -  dissipation 

Of  these  four  terms,  the  convection  (pu(3k/3x)  +  pv(3k/3y)J  and  production  [r(3u/3y)J 
terms  are  identically  those  derived  from  the  Navier-Stokes  equations  for  a  turbulent  flow, 
under  the  boundary-layer  assumptions,  along  with  the  assumptions  that  the  normal  stresses 
are  negligible  and  that  the  only  significant  tangential  stress  is  the  term  r  =  pu'v'.  The 
diffusion  and  dissipation  terms  are,  on  the  other  hand,  models  used  to  eliminate 
higher-order  correlations  of  fluctuating  terms  (about  which  little  is  known)  from  the 
equations. 

2.3.2  Dissipation 

There  seems  to  be  little  controversy  about  the  dissipation  term,  which  is  written 
as 

dissipation  =  a2pk3/,2/fk  (14) 

based  on  a  suggestion  of  Kolmogorov  (Ref.  17).  As  Reynolds  (Ref.  18)  points  out,  this 
term  does  not  represent  the  complete  dissipation  tensor  except  in  isotropic  flow,  and 
thus  should  be  referred  to  as  "isotropic  dissipation"  perhaps;  on  the  other  hand,  the  use 
of  the  term  to  describe  TKE  dissipation  in  general  is  nearly  universal. 

The  parameter  a2  in  Eq.  (14)  can  be  looked  upon  as  a  factor  usable  in  adjusting 
the  "isotropic  dissipation"  to  apply  to  nonisotropic  flows.  In  this  work  an  empirical 
formulation  for  this  parameter,  which  is  based  on  the  empirical  formulation  developed 
by  Peters  and  Phares  (Ref.  12),  has  been  used.  The  a2  relationships  used  in  this  work 
are  not  identical,  however,  to  those  used  by  Peters  and  Phares.  This  is  because  of  differences 
between  the  TKE  profile  development  predicted  by  the  method  and  the  profiles  assumed 
in  the  integral  method  of  Ref.  12. 


2.3.3  The  Turbulent  Reynolds  Number  and  Length  Scale 


To  describe  the  variation  of  a2,  a  "turbulent  Reynolds  number"  defined  by 


(15) 
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is  used  as  the  independent  variable.  This  parameter  is  a  function  of  x  only.  The  term 
em  represents  the  eddy  viscosity  at  the  maximum  shear  point  in  a  lateral  profile,  i.e., 

=  (VPmV(<5u/<Jy)m  (16) 

and  is  a  characteristic  length  scale,  effectively  the  lateral  width  of  the  turbulent  flow. 
This  scale,  which  is  also  used  in  the  dissipation  expression  (Eq.  (14)).  is  determined  by 
the  following  set  of  relationships: 

For  the  two-dimensional  shear  layer,  or  first  regime  of  a  jet  (sec  Fig.  1) 

£k  =  Ay  (17) 

where  Ay  is  the  lateral  distance  between  the  points  at  which  (u-ue)/(Uj-uc)  =  0.95  and 
(u-ue)/(uj-uc)  =  0.05,  for  profiles  which  are  not  fully  developed.  For  fully-developed 
profiles,  which  can  be  approximated  by  a  cosine  function, 

ek  =  1.57  Au/(«9u'<?v)m  (18) 

where  the  subscript  m  again  refers  to  the  position  in  a  lateral  profile  at  which  the  shear 
stress  is  a  maximum.  Equation  (18)  provides  a  characteristic  length  scale  determination 
which  is  compatible  with  the  cosine  function  velocity  profiles  used  in  Ref.  12,  and  which 
avoids  the  problem  of  defining  the  edge  of  an  asymptotic  profile.  Because  the  cosine 
function  is  not  always  a  good  approximation  to  the  true  computed  velocity  profile,  in 
practice  Eq.  (18)  is  only  used  in  place  of  Eq.  (17)  wherever  the  value  of  2*  determined 
by  it  satisfies  the  inequality  Ay  <  1.57  Ay.  For  a  two-dimensional  wake,  the  value 

of  Ay  to  be  used  in  Eq.  (17)  is  defined  as  the  lateral  distance  betwen  the  centerline 
and  the  point  at  which  (u-ucj/fupu,;)  =  0.05. 

For  the  second  regime  of  jets  (see  Fig.  1), 

Ek  =  2rH  09) 

where  rw  is  the  value  of  r  at  which  (u-uc)/(uc-uc)  =  0.5. 

In  the  definition  of  Rx,  Au  =  uj-ue  in  the  first  regime  of  jets,  or  the  two-dimensional 
shear  layer,  and  Au  =  uc-uc  in  the  second  regime  of  jets  and  in  wakes. 

The  variation  of  aj  with  Rx  was  determined  by  a  trial-and-error  technique  using 
the  variation  of  the  two-dimensional,  shear-layer  spreading  parameter  a  with  Mach  number 
(Test  Case  2  of  Ref.  1 1)  as  the  criterion.  Because  some  controversy  over  the  appropriate 
variation  of  a  with  M  in  this  case  exists,  it  should  be  pointed  out  that  for  the  most 
part  the  recommended  data  from  Ref.  1 1  were  used.  These  data  are  consistent  with 
the  observed  variation  of  velocity  potential  core  length  in  the  submerged  jet  as  a  function 
of  Mach  number  used  as  a  criterion  for  the  variation  of  a2  with  Rx  by  Peters  and  Phares 
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(Ref.  12),  The  expressions  which  resulted  from  the  trial-and-error  process  undertaken  in 
this  work  are: 


For  0  <  R-|-  <  185,  ilo  =  1-69 

For  18.i  <  Ft-|  360,  a,  -  0.46  +  0.00762  R  r 


For  Rt  >  360,  a,  =  3.20 

2.3.4  Modification  for  the  Effects  of  Density  Variation 


(20) 

(21) 

(22) 


Equations  (20)  through  (22)  describe  the  axial  variation  of  a; .  However,  further 
correction  is  necessary  (Ref.  12)  to  allow  lor  the  effects  of  density  gradients.  This  is 
obtained  by  multiplying  the  parameter  a?  obtained  from  Eqs.  (20)  through  (22)  by  the 
factor  (I/C,),  where,  for 


and  for 


Pr  | 

—  >  1 

c  *  =  0.984  +  0.0)6  - 

(23) 

p  jl 

1  Pjl 

ia< , 

Pit 

t:l  =  0-93  +  0.05ReTTe/RjT.rj 

(24) 

In  Eqs.  (23)  and  (24).  the  subscript  I  refers  to  conditions  at  the  origin  of  mixing, 
Rc  is  the  gas  constant  appropriate  to  the  outer  How,  and  Tje  is  the  total  temperature  of 
that  How.  Equation  (24)  is  used  when  the  jet  density  is  greater  than  the  tree-stream  density 
in  order  to  preserve  the  aj-R-p  relationship  for  the  supersonic  shear-layer  case,  for  which 
the  "jet"  side  is  taken  to  be  the  high-speed  stream.  Thus,  for  the  case  of  compressible 
but  identical-gas  mixing  problems,  the  value  of  C|  obtained  through  use  of  Eq.  (24)  is 
unity,  while  for  gases  with  different  molecular  weights  or  stagnation  temperatures  a 
nonunity  cj  is  obtained. 

2.3.5  Significance  of  the  Dissipation  Term 

There  are  three  points  at  which  empirical  information  enters  the  modeling  of  the 
turbulent  kinetic  energy  equation:  in  the  diffusion  term,  in  the  dissipation  term,  and  in 
the  relationship  between  the  turbulent  shear  stress  and  the  turbulent  kinetic  energy.  In 
this  paper,  the  most  extensive  empiricism  has  been  applied  to  the  dissipation  term,  rather 
than  to  the  diffusion  term  or  the  r-k  relationship,  because,  as  Bradshaw  and  others  (Ref. 
15)  point  out  with  reference  to  the  two-dimensional  boundary  layer,  the  dissipation 
parameter  is  distinctly  the  most  critical  term  in  the  modeling.  One  of  the  major  advantages 
of  a  simple  TKE  analysis  such  as  that  deicribed  in  this  report  is  that  the  facility  with 
the  appropriate  empiricism  may  be  introduced.  Thus,  while  the  more  rigorous  analyses 
such  as  those  proposed  by  Donaldson  and  Rosenbaum  (Ref.  10),  or  with  somewhat  greater 
simplification  by  llanjalic  and  Launder  (Ref.  19).  may  indeed  be  more  appropriate  than 
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the  present  analysis  for  situations  in  which  rapid  changes  in  the  turbulent  structure  occur, 
they  are  considerably  more  difficult  to  use  for  predicting  the  general  character  of  a  variety 
of  fiee  turbulent  Hows.  This  situation  is  caused  by  the  complicated  interrelationships  that 
exist  between  the  various  parameters  that  still  must  be  modeled  in  these  more  rigorous 
analyses,  as  well  as  the  remoteness  of  these  parameters  from  readily  available  experimental 
information.  The  present  analysis  is  intended  to  provide  general  predictive  (as  well  as 
correlative)  capability,  and  thus  some  approximation  of  the  details  of  the  turbulent 
structure  is  accepted  in  order  to  provide  a  more  general  applicability  of  the  mean  flow 
field  prediction. 

2.3.6  Diffusion  of  Turbulent  Kinetic  Energy 

The  remaining  term  to  be  modeled  in  the  TKF.  transport  equation  is  the  diffusion 
term,  which  is  here  approximated  by  a  gradient  diffusion  model  (Eq.  (6»  in  order  to 
retain  the  formal  similarity  of  the  equations,  discussed  in  Section  2.1.  There  is  some 
controversy  in  the  literature  over  this  point,  since  it  is  possible  to  avoid  a  gradient  model, 
as  was  done  by  Bradshaw  and  others  (Ref.  15)  for  the  two-dimensional  boundary  layer, 
and  by  Laster  (Ref.  20)  for  planar  and  axisymmotric  free-turbulent  flows.  The  significance 
of  the  approach  used  by  Bradshaw  and  others  and  by  Laster  is  that  the  TKE  equation 
becomes  hyperbolic,  which  both  allows  the  use  of  the  method  of  characteristics  for  the 
numerical  solution  of  the  problem  and.  more  significantly,  implies  that  turbulence  effects 
themselves  propagate  along  characteristic  lines  in  a  flow.  While  there  is  some  evidence 
that  this  type  of  behavior  is  experimentally  found  in  some  turbulent  flows  (see  Reynolds 
(Ref.  18)  for  a  discussion),  experimental  attempts  to  determine  the  nature  of  the  diffusion 
process  in  general  have  not  resulted  in  any  conclusive  evidence  in  favor  of  either  hypothesis. 
In  the  present  work,  the  TKE  is  assumed  to  diffuse  in  the  same  manner  as  any  other 
scalar  in  the  flow  field,  as  is  assumed  (albeit  with  a  different  model  for  the  r-k  relation) 
by  Rodi  and  Spalding  (Ref.  9).  This  has  the  advantage  of  retaining  the  same  form  for 
all  of  the  governing  equations. 

2.3.7  Turbulent  Prandtl  Number  for  Kinetic  Energy 

A  gradient  diffusion  model  requires  the  specification  of  a  diffusion  coefficient.  In 
common  with  the  other  scalar  transport  equations  the  TKE  diffusion  coefficient  is  taken 
to  be  related  to  the  momentum  diffusion  coefficient  e:  the  TKE  "Prandtl  number"  which 
provides  the  constant  of  proportionality  is  taken  to  be  constant  throughout  the  flow  field, 
with  a  value  of  0.7.* 

2.4  SUMMARY 

This  analysis  of  the  general  free-turbulent  mixing  problem  involves  the  simultaneous 
solution  of  the  turbulent  transport  equations  (Eqs.  (1)  through  (6)).  In  order  to  solve 
these  equations,  the  shear  stress  model  given  by  Eq.  (9)  is  used,  in  which  a  parameter 


♦The  solutions  have  been  found  to  be  relatively  insensitive  to  realistic  variation  in  the  Prandtl  number  (Ref.  34). 
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aj  appears.  This  parameter  is-  obtained  in  different  regions  of  the  flow  through  the 
relationships  given  by  Eq.  ( 1 0)-(  1 3).  To  obtain  the  turbulent  kinetic  energy  (TKE)  necessary 
in  Eq.  (9),  the  TKE  transport  equation  (Eq.  (6)).  is  solved.  Equation  (6),  however,  involves 
models  for  the  TKE  diffusion  and  dissipation.  Of  these,  the  latter  term  is  the  more 
important,  and  a  general  formulation  for  it  involves  the  definition  of  a  turbulent  Reynolds 
number,  Eqi  (15),  a  turbulent  length  scale,  Eqs.  (17)  through  (19),  and  a  dissipation 
coefficient.  a2,  Eqs.  (20)  through  (24). 

With  all  the  necessary  algebraic  relationships  defined,  the  system  of  equations  is  closed, 
and  a  numerical  solution  can  be  carried  out  subject  to  the  specification  of  the  initial 
and  boundary  conditions.  In  the  following  sections  an  outline  of  the  numerical  procedure 
will  be  given,  followed  by  a  discussion  of  the  specification  of  the  appropriate  initial  and 
boundary  conditions.  A  demonstration  of  the  abilities  of  the  method  will  then  be  exhibited, 
followed  by  a  user’s  manual  to  the  program  and  a  FORTRAN  listing. 

SECTION  III 

THE  NUMERICAL  PROCEDURE 
3.1  GENERAL  DESCRIPTION 

The  numerical  framework  used  in  the  program  delineated  in  this  report  is  a 
modification  of  the  procedure  devised  by  Patankar  and  described  in  Ref.  13.  In  this 
technique  the  governing  equations  are  transformed  into  nondimensional  stream  function 
coordinates,  allowing  the  lateral  grid  size  to  be  automatically  adjusted  to  account  for  the 
growth  of  the  shear  region,  which  can  be  quite  rapid  in  a  free-mixing  flow.  This  grid 
adjustment  feature  and  the  successive  substitution  method  used  to  evaluate  the  solution 
matrix  at  each  downstream  step  combine  to  produce  an  efficient,  fast-running  program. 

Two  major  modifications  were  made  to  the  Patankar  technique  to  produce  the 
program  described  in  this  report.  Relatively  straightforward,  the  first  is  the  incorporation 
of  the  turbulent  kinetic  energy  equation  and  the  shear  stress  model  used  in  this  work. 
The  second  major  modification  was  made  to  ensure  the  satisfaction  of  the  appropriate 
integral  conservation  equation,  which  requires,  in  the  case  of  a  zero  pressure  gradient  free 
mixing  flow,  that  the  momentum  excess  (or  defect)  integral  be  a  constant  everywhere 
in  the  flow.  The  original  technique  did  not  satisfy  this  requirement,  for  reasons  which 
are  described  in  Part  II  of  Ref.  13. 

In  the  remainder  of  this  section  a  brief  description  of  the  development  of  the 
finite-difference  procedure  will  be  given.  Because  the  procedure  is  described  at  length  by 
Patankar,  only  the  highlights  of  the  method  will  be  discussed  here;  it  is  recommended 
that  the  user  of  this  program  become  familiar  with  the  contents  of  Ref.  13. 
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3.2  COORDINATE  TRANSFORMATION 

To  put  Eqs.  (1)  through  (6)  into  a  form  suitable  for  rapid  numerical  computation, 
a  stream  function  ip  is  defined,  which  satisfies  the  continuity  Eq.  (1)  identically.  Thus, 
at  constant  x 


purdy  =  d 0  (25) 

where  r  =  rj  +  y  cos  Q  (Fig.  3).  Next,  a  nondimensional  stream  function,  co,  is  defined, 
such  that 


(26) 


where  E  and  I  represent  the  external  and  inner  boundaries  of  the  mixing  zone,  respectively 
(see  Fig.  3).  Note  that  from  Eq.  (26),  0  <,  co  <  1.  The  term  (i//e  -  ^i)  represents  the 
mass  flow  in  the  mixing  region  at  a  given  x,  i.e., 

2'e 

^e:  ~  't'l  =  J  Purd y  (27) 


After  conversion  to  x-oj  coordinates,  each  of  the  Eqs.  (2)  through  (6)  can  be  written 
in  the  standard  form 


§*  +  <a+b*>)f2 

0X  0(0 


(28) 
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where 


(29) 

(fj-mg-  fjiiip/ (^E  - 

(30) 

y2a  -  ^,)2 

(31) 

In  these  expressions,  m"i  =  p\\\  represents  the  mass  flux  per  unit  area  entrained  into 
the  mixing  region  through  the  I  boundary,  while  riV't  =  -Peve  is  the  same  quantity 
evaluated  at  the  E  boundary.  The  definitions  of  00  and  d  depend  on  the  dependent  variable 
and  are  given  in  Table  II. 


TABLE  II 

SOURCE  TERMS  IN  THE  TRANSFORMED  EQUATIONS 


3.3  SPECIFICATION  OF  THE  AUXILIARY  FUNCTIONS 
3.2.1  Entrainment 

It  is  in  the  evaluation  of  the  terms  (^e  -  '/h),  tbitT'e,  and  rim"i  that  substantial 
differences  have  evolved  between  the  formulations  used  by  Patankar  and  Spalding  (Ref. 
13)  and  those  used  in  this  work.  In  the  case  of  the  entrainment  rates,  the  differences 
have  resulted  from  the  use  of  a  shear  stress  model  which  does  not  incorporate  a  mixing 
length  concept,  while  the  difference  in  the  evaluation  of  (i^e  -  tf'i)  results  from  the 
requirement  that  the  integral  momentum  equation  be  exactly  satisfied. 
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The  terms  rEm"E  and  rim"i  are  evaluated  from  the  axial  momentum  equation,  which 
can  be  written 


d  y2gpuf  da 

dci  doi_ 
da/dco 


At  the  I  and  E  boundaries,  9u/9<d  «  0  and  Eq.  (32)  cannot  be  used  to  obtain  the 
entrainment  rates.  However,  this  equation  may  t>e  applied  just  inside  the  edges  of  the 
mixing  region,  for  example,  at  cj  =  0.05  and  at  gj  =  0.95.  Consider  the  evaluation  at 
the  outer  edge,  where  cj  =  cjb  =  0.95.  The  only  term  in  Eq.  (32)  which  cannot  be  directly 
evaluated  (assuming  a  known  axial  pressure  gradient  is  duB/dx,  since  at  this  point  the 
downstream  velocity  profile  is  not  known.  However,  if  it  is  required  that  at  the  downstream 
step  the  velocity  be  some  desired  value,  for  example,  that  Ub  =  U|  +  0.99  (ue  -  uj), 
then 


l“‘l 


r fiHi  +  *  ^lmI^  = 


^E-^I)/duB  1 

“  -7T-m  '■  I  ~T—  +  - 


(da/da)B  \  d*  Pbub 


$ 


(32) 


dx 


(33) 


where  uB  is  the  value  of  velocity  acutally  achieved  at  co  =  coB  in  the  presently  known 
profile.  A  similar  prescription  can  be  formulated  for  the  inner  edge,  if  necessary,  and 
the  entrainment  rates  evaluated  from  Eq.  (32).  If  one  boundary  is  a  wall  or  an  axis  of 
symmetry,  Eq.  (32)  is  evaluated  at  the  opposite  boundary;  if  both  boundaries  are  free, 
the  entrainment  rates  must  be  obtained  at  both  boundaries  simultaneously. 


In  practice,  the  most  important  effect  of  the  entrainment  rate  specification  is  on 
the  size  of  the  lateral  grid  spacing.  Clearly,  if  the  grid  spacing  becomes  too  large,  accuracy 
is  lost,  because  too  many  of  the  computational  points  will  be  outside  the  mixing  region 
of  interest.  Conversely,  if  the  spacing  becomes  too  small,  there  will  not  be  enough  points 
available  to  completely  cover  the  region  of  interest.  In  order  to  keep  the  grid  spacing 
small  enough  for  accuracy,  many  controls  can  be  applied:  perhaps  the  most  effective  is 
to  monitor  the  nondimensional  velocity  profile  <£u  =  (u-uc)/(uc-ue)  at  the  edges.  Then, 
if  at  the  preceding  station  this  value  has  decreased  below  some  previously  defined  ^u*, 
the  current  entrainment  rate  may  be  multiplied  by  <f,u/<&u*>  Increasing  an  emrainment 
rate  which  is  too  low  is  easily  accomplished  by  alteration  of  the  value  of  cjb  at  which 
the  entrainment  is  computed. 


Too  much  emphasis  should  not  be  placed  on  the  method  of  specification 'of  the 
entrainment.  It  may  be  obtained  from  any  of  the  transport  equations  or  in  any  other 
way  consistent  with  the  physics  of  the  problem,  since  an  incorrect  specification  is  always 
clearly  apparent  from  the  velocity  profiles  obtained  from  the  solution.  The  techniques 
used  in  the  program  listing  included  in  this  report,  in  conjunction  with  the  profile  shape 
controls,  allow  a  quite  general  determination  of  the  entrainment  rates;  in  fact,  all  of  the 
computations  to  be  described  in  the  following  section  have  been  carried  out  using  the 
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same  specification  technique.  Should  the  specification  fail,  however,  any  alternate  technique 
which  allows  the  computation  to  proceed  smoothly  is  acceptable. 


3,3.2  The  Mixing  Region  Mass  Flux 


In  Patankar's  and  Spalding's  original  program  (Ref.  13)  the  mass  flux  (0E  -  i/q)  at 
the  downstream  station  is  evaluated  as  the  sum  of  -  0q)  at  the  upstream  station 
plus  the  entrained  mass  flux  over  the  step.  However,  inaccuracies  in  the  finite-difference 
formulation  at  the  edge  regions  (which  are  discussed  in  Part  II  of  Ref.  13)  introduce 
a  small  error  in  the  determination  of  the  coefficients  in  the  finite-difference  version  of 
the  equations  of  motion,  which  is  reflected  in  a  cumulative  inaccuracy  in  the  momentum 
integral  at  each  station.  Referring  to  Fig.  1,  the  momentum  excess  (or  deficit)  integral 
at  any  station  can  be  written 


rr  1 E 

PjUj(Uj  -  ue)  —  +  J  pu(u-ue)rdy  =  constant 


(34) 


which,  under  the  transformation  of  coordinates  described  above,  becomes 

r?  1 

Pjuj(uj  -  ue)  —  +  (^E  -  ^[)  f  (u  -  ue)d(u  »  constant  (35) 

Thus,  Eq.  (35)  can  be  used  to  evaluate  (0E  -  0|)  subject  to  the  condition,  valid  for 
a  constant-pressure  free  turbulent  flow,  that  Eq.  (35)  is  satisfied.  Note  that  Eq.  (35)  is 
not  appropriate  in  the  case  of  nonzero  axial  pressure  gradient;  however,  the  appropriate 
modifications  may  easily  be  introduced  provided  that  dp/dx  is  known  externally  to  the 
mixing  calculations. 

3.4  DEVELOPMENT  OF  THE  DIFFERENCE  EQUATIONS 

The  details  of  the  numerical  procedure  developed  by  Patankar  are  described  in  Ref. 
13;  it  is  an  implicit  procedure  similar  to  the  Crank-Nicholson  technique.  However,  Patankar 
and  Spalding's  formulation  of  the  difference  equations  varies  sufficiently  from  other 
formulations  to  warrant  a  brief  description. 

Consider  the  finite-difference  grid  shown  in  Fig.  4.  The  points  U-,  U,  and  U+ 
correspond  to  the  three  values  up.,  cup ,  and  gju+;  the  corresponding  downstream  points, 
at  equal  values  of  ou,  are  denoted  by  D-,  D,  and  D+.  It  is  assumed  that  between  grid 
points  the  dependent  variables  0  vary  linearly  with  co,  and  that  along  the  x-coordinate 
between  xE-  and  xp  the  value  of  0  is  0u  except  at  x  =  xp,  where  it  becomes,  with 
a  step  change,  0D.  The  values  of  0  at  x  =  Xp  are  all  known,  and  the  values  of  <f>  at 
x  =  xp  are  to  be  solved  for  simultaneously,  using  these  values.  If  a  linear  variation  in 
0  between  xp  and  xp  were  assumed,  the  method  would  correspond  to  the  Crank-Nicholson 
technique. 
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Fig.  4  Finite-Difference  Grid 


To  formulate  the  difference  equation,  the  values  of  the  derivatives  are  obtained  as 
mean  values  integrated  over  the  control  volume  shown  crosshatched  in  Fig.  4.  The  points 
UU-,  UU+,  DD-,  and  DD+  are  all  evaluated  halfway  between  their  respective  node  points 
on  an  x-to  grid.  In  keeping  with  the  formulation  of  the  difference  equation  as  a  "miniature 
integral  equation,"  the  fundamental  derivatives  are  expressed  as 


36 

<5x 


do>  dx 


1 - JDD+ 

w  D  D+  —  W  D  D—  W  D  u_  doj 


(36) 


The  integrals  are  evaluated  using  the  assumed  linear  profiles  between  grid  points; 


thus 


Xn  (J  n  n  W  n  n _  ^  » T  _ 


“dd+  /d 


j*  D  D+ 


dx  — -  -j  j  d% 
WDD-  XU 


DD- 


0>j 


^rDD  uut 

=  J  —  <5|-)dct>  +  J 

“dd-  wdd  (37) 


DD+ 


DD 


In  general,  for 


wDD-  -  w  S  WD  ^  ~  ^DD-  +  ^  |(w  — 
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and  for 


^  ^  &>QQ+  0  “  +  K  ry{0J  —  CO  |}) 


where 


~  (<£d-“  ^DD-^  -  bDdJ'  etc‘ 

Substituting  these  expressions  into  Eq.  (37)  leads  to  the  expression 


d<f)  ~  ^D-"^U^D"aD-*  3  ^D+ ~  <^ll+^ea13+ 

d\  4(xd  -  XU)(£JD+  -  <UD_)  4  (xD  -  Xy)  4(xd  —  xI()(&)D+  -  «aD_) 


(38) 


Similarly,  for  90/ 9a>,  evaluation  of  Eq.  (36)  yields 


B±  %  l  ,“»>*■ 

ddi  aDD+_(aDD-  cjdd_ 


d  (f> 


so  that 


90  ^  Pp+  ^p- 
do>  wd+-"d- 


SECTION  IV 

A  DEMONSTRATION  OF  THE  METHOD 
4.1  SPECIFICATION  OF  INITIAL  AND  BOUNDARY  CONDITIONS 


(39) 


The  analysis  outlined  in  the  preceding  sections  allows  the  rapid  computation  of  a 
variety  of  free  turbulent  mixing  problems  which  are  themselves  differentiated  by  the  initial 
and  boundary  conditions  involved.  The  initial  conditions  which  need  to  be  specified  involve 
beginning  profiles  of  the  velocity,  concentration,  total  enthalpy,  and  turbulent  kinetic 
energy,  with  specification  of  the  concentration  profile  necessary  only  in  a  two-gas  problem 
and  specification  of  the  total  enthalpy  profile  necessary  only  in  nonisoenergetic  flow.  The 
appropriate  boundary  conditions  depend  on  the  nature  of  the  problem.  Thus,  at  an  axis 
of  symmetry,  all  lateral  gradients  are  assumed  to  be  zero,  while  at  a  free  boundary,  uniform 
values  of  velocity,  concentration,  total  enthalpy,  and  turbulent  kinetic  energy  are  specified. 
No  fundamental  changes  in  the  program  are  necessary  to  accommodate  prescribed  axial 
variation  of  each  of  these  parameters. 

Because  of  the  coordinate  system  used  in  this  analysis,  the  initial  profiles  involved 
cover  only  the  initial  viscous  region.  In  general,  if  a  calculation  is  started  from  the  nozzle 
exit  plane,  power-law  initial  velocity  profiles  are  used,  i.e..  (u/u B)  =  (y/S),/ln.  where  ub 
represents  either  the  jet' or  outer  stream  velocity.  These  power-law  profiles  are  patched 
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together  at  the  nozzle  lip  assuming  a  nonzero  but  small  velocity,  typically  0.01  Uj  (see 
Fig.  5).  Uniform  (top-hat)  profiles  of  C  and  H  are  generally  assumed,  if  required. 


a.  Velocity  b.  Enthalpy  or  c.  Eddy  Viscosity 

Concentration 


Fig.  5  Initial  Condition  Definition  Sketch 
4.1.1  Specification  of  the  Initial  Kinetic  Energy  Profile 

Because  the  turbulent  kinetic  energy  equation  is  used  in  this  analysis  to  obtain  the 
turbulent  shear  stress  in  the  flow  field,  specification  of  the  initial  turbulent  kinetic  energy 
profile  is  required.  While  the  solution  downstream  is  not  overly  sensitive  to  this  initial 
specification,  some  care  has  to  be  exercised  in  providing  it.  and  since  initial  turbulent 
kinetic  energy  levels  are  not  ordinarily  part  of  the  specification  of  a  given  engineering 
problem,  a  general  means  of  developing  approximate  initial  TKE  profiles  has  been 
developed. 

Two  techniques  for  obtaining  initial  TKE  profiles  have  been  used  with  considerable 
success  in  the  demonstration  of  the  method  which  will  follow.  If  the  calculation  is  to 
begin  at  the  nozzle  lip.  the  eddy-viscosity  profiles  developed  for  compressible  turbulent 
boundary  layers  by  Maise  and  McDonald  (Ref.  21)  are  used  to  obtain  the  shear  stress 
and  thus  the  turbulent  kinetic  energy  through  use  of  Eq.  (9).  Alternatively,  if  the 
computation  is  to  begin  downstream  of  the  nozzle  lip.  the  Prandtl  eddy  viscosity  model 
(Ref.  4)  is  used  to  obtain  the  initial  shear  stress  and  hence  the  kinetic  energy  level.  The 
constants  used  for  this  model  are  shown  in  Table  111.  (This  level  may  also  be  obtained 
in  some  cases  from  experimental  data.)  One  or  the  other  of  these  techniques  has  been 
used  for  all  of  the  calculations  to  be  presented  in  the  remainder  of  this  section;  the  specific 
model  used  in  each  case  is  described  in  Table  III. 
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(Ref.  11) 

Description 

Start 

Initial  Conditions 

No. 

Point  (CM) 

u 

c 

H 

k 

i 

Incompressible  2-D  shear 
layer 

x  =  0 

a 

M/A 

N/A 

b 

2 

Compressible  2-D  shear 
layer 

x  •  0 

a 

N/A 

c 

b 

3 

Variable  density  2-D  shear 
layer 

x  =  0 

a 

c 

c 

b 

4 

Incompressible  2-D  shear 
layer  (Lee) 

x  =  0 

d 

N/A 

N/A 

b 

5 

Compressible  2-D  shear 
layer  (Hill  and  Page) 

K 

II 

to 

d 

N/A 

c 

e 

6 

Circular  jet  (Maestrello 
and  McDaid) 

x/D  ■  1 

d 

N/A 

N/A 

f 

7 

Supersonic  circular  jet 
{Eggers) 

x  ®  0 

d 

N/A 

c 

b 

8 

Compressible  circular  jet 
(Heck) 

x/D  =  2.B 

d 

N/A 

d 

f 

,  S 

Coaxial  air  jets  with 

He  trace  (ForatalL) 

x/D  *  0 

d 

c 

c 

b 

10 

Coaxial  H2-alr  jets 
(Chriss) 

x/D  *  0 

d 

d 

d 

t 

11 

Compressible  coaxial 
air- air  jets  (Eggers  and 
Torrence) 

x/D  '  0 

d 

N/A 

c 

b 

12 

Compressible  coaxial 

H2-air  jets  (Eggers) 

x/D  «  0 

d 

c 

c 

b 

13 

Plane  2 -stream  jets 
(Bradbury) 

x  3  0 

a 

N/A 

N/A 

b 

M 

Incompressible  axl sym¬ 
metric  wake  (Chevray 
and  Kovasznay) 

x/D  ■  0 

d 

N/A 

N/A 

b 

15 

Incompressible  axisym- 
metric  wake  (Chevray) 

x/D  •  0 

d 

N/A 

N/A 

f 

16 

Compressible  plane  wake 
(Demetriades) 

x  3  0.91 

d 

N/A 

c 

K 

17 

Compressible  axisym- 
metric  wake  (Demetriades) 

x  =  6.  74 

d 

N/A 

c 

h 

IB 

Circular  jet  (Wyynanaki 
and  Fiedler) 

x/D  •  1 

i 

N/A 

N/A 

f 

19 

Compressible  circular  ]et 
(Heck) 

x/D  •  2.  8 

d 

N/A 

c 

f 

20 

Coaxial  air-air  Jet  (Paulk) 

x/D  >2.3 

d 

a 

a 

f 

21 

Coaxial  Hj-air  jets 
(Chriss) 

x/D  *  2.6 

d 

a 

a 

( 

22 

Coaxial  compressible 

H2~air  jets  (Eggera) 

x/D  *  0 

d 

_ 

c 

c 

b 

Legend: 


a.  Power  >law  profile,  n  =■  1/7. 

b.  Maiae  and  McDonald  (Ref.  21)  eddy  viscosity  profile  and  level. 

c.  Constant  value  baaed  on  data  sheet  (Ref.  11). 

d.  Experimental  values  from  data  sheet  (Ref.  11). 

e.  Maiae  and  McDonald  (Ref.  21)  eddy  viscosity  on  high-velocity 
side,  constant  eddy  viscosity  (equai  to  one-half  adjacent 
Maine  and  McDonald  value)  for  low  speed  side.  (Actual  value 
unimportant,) ) 

f .  Prandtl  (Ref.  fe)  eddy  viscosity,  kp  c  0.005. 

g.  Constant  eddy  viscosity,  level  based  on  data  of  Ref.  22. 

h.  Constant  eddy  viscosity,  level  based  on  R*.  data  of  Ref.  23. 

1 .  Data  from  Ref.  24* 
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4.2  GENERAL  DESCRIPTION  OF  THE  DEMONSTRATION 

The  data  chosen  for  the  demonstration  of  the  method  is  that  used  in  the  1972 
NASA-Langley  Working  Conference  on  Free  Turbulent  Shear  Flows  (Ref.  11).  These  data 
were  chosen  for  two  reasons:  (1)  they  represent  a  compilation  of  a  variety  of  free  mixing 
experiments  against  which  a  variety  of  models  were  tested,  and  (2)  the  present  model 
represents  a  modification  and  improvement  of  a  model  presented  at  that  conference  (Ref. 
2),  and  thus  prediction  of  the  same  cases  will  allow  a  direct  comparison  of  the  present 
predictions  with  those  of  Ref.  2.  It  should  be  emphasized  that  the  program  listed  in 
Appendix  A  was  used  to  compute  all_of  the  flows  shown  here;  there  were  no  changes 
made  in  the  empirical  information  used  or  to  the  program  structure  in  order  to  handle 
any  of  the  computations. 

Four  types  of  flow  fields  were  considered  in  the  data  package  for  the  1972  NASA 
Conference:  the  two-dimensional  shear  layer,  circular  jets  into  still  surroundings,  two-stream 
jets  (both  axisymmetric  and  two-dimensional)  and  axisymmetric  and  two-dimensional 
wakes.  In  the  following  subsections  the  predictions  made  by  the  present  model  for  each 
of  these  flow  classes  will  be  described. 

4.3  THE  TWO-DIMENSIONAL  SHEAR  LAYER 

Cases  1-3  of  the  NASA  Conference  involved  the  prediction  of  the  behavior  of  the 
two-dimensional  shear  layer  spreading  rate  as  a  function  of  velocity  ratio,  density  ratio 
and  Mach  number.  The  spread  parameter  a  was  computed  from  the  standard  definition 
used  at  the  conference: 


<r  -  1.855  (Xj-XjV^-Yj)  (40) 

Here,  Y2  represents  the  lateral  distance  between  the  points  at  which  (u-uc)/(uj-ue) 
=  0.9  and  (u-ue)/(uj-ue)  =  0.1  at  X2,  with  Yi  representing  a  similar  distance  at  Xj.  The 
subscripts  j  and  e  represent  the  primary  and  secondary  streams. 

In  the  version  of  this  computational  approach  reported  at  the  NASA  Conference 
(Ref.  2),  a  special  profile  of  the  parameter  aj  was  used  which  was  applicable  to  only 
the  two-dimensional  shear  layer.  Subsequent  to  the  presentation  of  Ref.  2,  the  use  of 
a  special  ai  function  was  found  to  be  unnecessary.  In  the  present  work  the  same  ai 
function  (described  by  Eq.  (10)  through  (13)  has  been  used  for  all  computations. 

Figure  6  shows  that  the  prediction  of  the  variation  in  the  spreading  parameter  a 
with  velocity  ratio  obtained  with  the  present  model  closely  follows  the  classic  expression 


1-Vni 

1  +  ue/u. 


(41) 


which  is  widely  accepted  as  the  proper  relationship  for  growth  rate  as  a  function  of  velocity 
ratio  (Birch,  in  Ref.  11). 
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Seconder y/Prlaary  Steam  Velocity/Ratio,  ue^uj 

Fig.  6  Predicted  Effect  of  Velocity  Ratio  on  Two-Dimensional 
Shear- Layer  Spread  Rate,  NASA  Case  1 

The  variation  of  the  spreading  parameter  a  with  Mach  number  is  still  a  controversial 
issue,  although  the  proper  prediction  of  this  spread  parameter  is  crucial  to  the  correct 
prediction  of  the  length  of  the  potential  core  in  a  supersonic  jet.  However,  a  recommended 
set  of  data  for  the  variation  of  o  with  Mach  number  emerged  from  the  NASA  Conference, 
and  these  data,  together  with  the  prediction  made  by  the  present  model,  are  shown  in 
Fig.  7.  As  described  in  Section  2.3.3.  these  computations  were  used  to  develop  the  a2*Rj 
expressions,  Eqs.  (20)  through  (22).  Except  for  the  point  at  M  =  1,  the  agreement  is 
quite  good.  (There  is  some  evidence,  primarily  sonic  jet  core  lengths,  that  the  spread  rate 
data  point  at  M  =  1  represents  a  value  of  a  which  is  too  low;  thus,  in  developing  Eqs. 
(20)  through  (22),  which  in  effect  control  the  shape  of  the  curve  of  oo/o  versus  M,  more 
attention  was  paid  to  proper  prediction  of  the  data  at  M  =  2  and  higher.) 

The  predicted  variation  of  a  with  density  ratio  is  shown  on  Fig.  8.  Here  the  correct 
physical  behavior  is  again  controversial,  and  there  are  no  data  available  which  are  directly 
comparable  (i.e.,  with  uc/uj  =  0.2)  with  these  predictions.  However,  the  data  discussed 
in  Ref.  1 1  at  other  velocity  ratios  appear  to  show  a  smaller  variation  of  a  with  density 
ratio,  than  that  observed  with  Mach  number.  This  is  supported  by  these  calculations,  at 
least  for  1/8  <  pj/pc  ^  14.  Two  curves  illustrating  the  effect  of  Prandtl  number  change 
are  shown.  The  calculations  were  made  with  a  density  variation  created  through  a 
temperature  difference  between  the  streams.  Within  the  framework  of  the  present  analysis, 
and  assuming  a  Lewis  number  of  unity,  the  same  effects  would  be  observed  if  the  density 
difference  were  caused  by  molecular  weight  differences  between  the  streams. 
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Kach  Number  of  Primary  Stream 

Fig.  7  Predicted  Variation  of  Spread  Parameter  with  Mach  Number, 
NASA  Case  2 


Primary  Streaa  Density/Secondary  Stream  Density,  pj/pe 


Fig.  8  Variation  of  Spreading  Parameter,  with  Density  Ratio, 
NASA  Case  3 
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These  three  cases  so  far  discussed  have,all  involved  fully  developed  free  shear  layers. 
In  the  next  two  cases,  computations  have  been  made  of  the  developing  free  shear  layer. 
Figure  9  shows  a  comparison  of  the  predicted  development  of  an  incompressible  two-stream 
free  shear  layer  with  experimental  data,  while  Fig.  10  provides  the  same  comparison  for 
a  compressible  free  shear  layer.  The  agreement  in  both  cases  is  excellent.  In  common 
with  the  other  boundary  layer  solutions,  this  computational  technique  does  not  satisfy 
the  proper  boundary  conditions  for  these  nonsymmetric,  two-dimensional  flows.  The  effect 
of  this  is  that  the  spatial  orientation  of  the  computed  profiles  is  not  known,  and  so  for 
the  comparison  the  experimental  and  theoretical  profiles  have  been  matched  at  the  point 
at  which  the  velocity  is  equal  to  the  average  velocity  of  the  two  streams.  It  should  also 
be  noted  that  because  the  numerical  method  utilizes  a  variable  x-step,  the  axial  locations 
of  the  predicted  and  measured  profiles  will  not  agree  exactly  in  general. 


Fig.  9  Comparison  of  Experimental  and  Theoretical  Velocity  Profiles 
for  NASA  Case  4  (Incompressible  Two-Dimensional  Shear 
Layer,  u2/u1  =  0.357) 
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'  1  ' _ ■  ■  1 

0  0.2  0.4  at  0.8  1.0 

U/U] 


Fig.  10  Comparison  of  Theoretical  and  Experimental  Profiles,  Case  5, 
Compressible  (M  =  2)  Two-Dimensional  Shear  Layer 

4.4  THE  CIRCULAR  JET 

The  prediction  of  the  asymptotic  behavior  of  an  incompressible  circular  jet  (NASA 
Conference  Case  18)  is  shown  in  Fig.  11,  from  which  it  can  be  seen  that  although  the 
prediction  of  the  centerline  velocity  is  proportional  to  (x/D)'1 ,  as  is  required  by  similarity 
considerations,  the  proportionality  constant  does  not  agree  exactly  with  that  which  may 
be  obtained  from  the  data  of  Albertson,  et  al.  (Ref.  25).  As  discussed  in  Ref.  2,  the 
data  proposed  for  the  comparison  in  NASA  Case  18  do  not  fully  represent  an  asymptotic 
jet.  For  this  reason  the  data  of  Ref.  25  have  been  used  for  this  comparison.  As  can  be 
seen  from  Fig.  12,  the  width  scale  of  the  flow  is  also  too  small,  which  agrees  with  the 
overprediction  of  the  similar-region  centerline  velocity  illustrated  in  Fig.  11.  Comparing 
the  profiles  of  velocity  and  turbulent  kinetic  energy  with  the  data  used  for  Case  18  (from 
Wygnanski  and  Fiedler,  Ref.  26)  also  shows  that  the  predicted  turbulent  kinetic  energy 
is  too  low  compared  to  the  experiment  on  the  centerline.  This  would  indicate  that  the 
TKE  diffusion  is  somewhat  too  small  near  the  centerline  in  the  asymptotic  jet,  Other 
computations  of  tliis  flow,  not  shown  herein,  indicate  that  a  lateral  variation  of  the  TKE 
Prandtl  number  may  be  required  to  obtain  better  agreement. 
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Axial  Distance/Nozzle  Diameter,  x/D 

Fig.  11  Comparison  of  Prediction  with  Experimental  Data,  NASA  Case  18, 
Asymptotic  Circular  Jet 


Radlus/Axlal  Distance,  r/x 

Fig.  12  Comparison  of  Profile  Predictions  with  Experimental  Data, 

NASA  Case  18,  Asymptotic  Circular  Jet 

Circular  jet  predictions  were  also  required  for  NASA  Cases  6,  7,  8  and  19,  under 
varying  conditions.  The  performance  of  the  present  model  in  predicting  these  cases  is 
shown  in  Figs.  13-16.  The  jet  Mach  number  for  these  cases  ranged  from  0.6  to  2.22, 
and  as  these  figures  show,  the  present  model  provides  satisfactory  predictions  in  each 
case.  There  appears  to  be  a  general  trend  in  these  predictions  to  underpredict  the  centerline 
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velocity,  which  may  be  traced  in  Cases  6  and  7  to  a  slight  underprediction  of  the  velocity 
potential  core  length. 


Axial  Distance  x 
Nozzle  Radius  rQ 

Fig.  13  Comparison  of  Prediction  with  Experimental  Data  for  M  =  0.6 


Circular  Jet,  NASA  Case  6 


Axial  Distance  x 
Nozzle  Radius  rQ 

Fig.  14  Comparison  of  Prediction  with  Experiment,  M  =  2.22  Jet,  NASA  Case  7 
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Fig.  15  Comparison  of  Prediction  with  Experiment,  NASA  Case  8,  Heated  M  =  0.7  Circular  Jet 


Fig.  16  Comparison  of  Prediction  with  Experiment,  NASA  Case  19,  Mj  =  1.4  Circular  Jet 
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4.5  TWO-STREAM  JETS 

The  eight  cases  of  two-stream  jet  data  selected  for  the  NASA  Conference  can  be 
divided  into  three  subgroups:  momentum  transport  alone;  momentum  and  energy  transport; 
and  momentum,  energy,  and  mass  transport.  The  simplest  configurations  were  the 
two-stream,  two-dimensional  mixing  process  represented  by  Case  13,  and  the  coaxial 
two-stream  process  represented  by  Cases  9  and  20.  In  the  latter  two  cases  energy  and 
mass  tracers  were  included,  but  in  small  enough  quantity  that  the  overall  flow  field  was 
not  affected. 

As  Figs.  17  and  18  show,  the  prediction  of  both  the  two-dimensional  flow  of  Case 
13  and  the  coaxial  configuration  of  Case  9  is  excellent.  The  prediction  of  the  coaxial 
flow  of  Case  20  (Fig.  19)  is  not  good;  the  results  of  the  conference  showed  that  this 
flow  was  quite  difficult  for  most  of  the  models  presented  to  predict. 

The  supersonic  mixing  configuration  represented  by  NASA  Case  1 1  represented  a 
flow  with  momentum  and  energy  transfer.  This  particular  case  was  further  complicated 
by  the  presence  of  extremely  thick  boundary  layers  at  the  nozzle  exit,  and  the  wake-like 
character  of  the  flow  in  which  the  outer  stream  velocity  is  greater  than  that  of  the  jet 
velocity.  A  prediction  of  this  case  by  the  present  model  is  shown  in  Fig.  20,  which 
demonstrates  a  good  prediction  of  the  boundary-layer-dominated  mixing  process. 


Fig.  17  Comparison  of  Prediction  with  Experiment,  NASA  Case  13, 
Two-Dimensional  Jet  in  Moving  Stream,  uj/ue  =  3.29 
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Fig.  18  Comparison  of  Prediction  with  Experiment,  NASA  Case  9. 
Coaxial  Air  Jets 


Fig.  19  Comparison  of  Prediction  with  Experiment,  NASA  Case  20, 
Coaxial  Air  Jets 
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Fig.  20  Comparison  of  Prediction  with  Experiment,  NASA  Case  11, 
Axisymmetric  Jet  in  Moving  Stream,  Mj  =  0.90,  Me  =  1.30 


Mixing  with  momentum,  energy,  and  mass  transfer  was  represented  by  four  cases, 
and  predictions  of  the  present  model  for  all  of  these  cases  are  shown  in  Figs.  21  through 
24.  Cases  10  and  21  (Figs.  21  and  22)  involve  the  prediction  of  the  subsonic  mixing 
of  coaxial  jets  of  hydrogen  and  air.  Both  of  these  cases  involve  start  points  downstream 
of  the  nozzle  exit.  Because  the  start  region  in  both  cases  was  in  the  region  of  transition 
between  first  and  second  regime  behavior,  a  special  start  technique  was  used  for  them. 
This  special  technique  involved  the  use  of  Eqs.  (12)  and  (13)  in  the  development  of  the 
initial  kinetic  energy  profile  from  the  initial  eddy  viscosity  profile  used  as  a  start  condition, 
and  the  subsequent  use  of  Eqs.  (12)  and  (13)  to  obtain  the  shear  stress  from  the  TKE 
profile  even  where  Uc  >  0.9  Up.  The  reason  for  this  alteration  of  the  standard  procedure 
is  that  the  late  start  point  (in  terms  of  flow  field  development)  did  not  allow  the  buildup 
of  centerline  turbulent  kinetic  energy  to  occur,  The  calculation  normally  predicts  such 
a  buildup  in  the  transition  region.  The  usejof  Eqs.  (12)  and  (13)  to  increase  the  TKE 
level  near  the  centerline  over  that  which  the  assumed  constant  eddy  viscosity  would 
produce,  in  conjunction  with  Eq.  (11),  compensated  for  the  lack  of  centerline  buildup, 
at  least  in  a  crude  fashion. 

Figures  23  and  24  illustrate  the  predictions  of  the  model  for  NASA  Cases  12  and 
22.  Both  of  these  cases  were  dominated  by  thick  boundary  layers  and  a  strongly  wake-like 
character  produced  by  the  low  momentum  flux  of  the  central  hydrogen  stream.  In  both 
cases  the  initial  decay  rate  is  overpredicted,  but  the  qualitative  agreement  with  the  data 
for  Case  22,  in  which  the  centerline  velocity  decreases  below  the  outer  stream  velocity 
and  then  slowly  increases,  is  noteworthy. 
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Fig.  21  Comparison  of  Predictions  with  Experiment,  NASA  Case  10, 
Hydrogen-Air  Jets 
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Fig.  22  Comparison  Theory  with  Experiment,  NASA  Case  21, 
Hydrogen-Air  Jets 
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at  Nozzle  Exit 


lomparison  of  Predictions  with  Experimental  Data,  NASA  12,  Hydrogen-Air  Jets 


Axial  Distance /Nozzle  Disaster,  x/D 


Fig.  24  Comparison  of  Predictions  with  Experiment,  NASA  Case  22, 
Supersonic  Hydrogen-Air  Jets 
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4.6  WAKES 

Four  wake  flows  were  included  in  the  test  cases  of  Ref.  1 1 ,  including  two-dimensional 
and  axisymmetriCj  incompressible  and  compressible  flows.  The  predictions  of  the  present 
model  for  these  cases  are  shown  in  Figs.  25  through  30. 


Axial  Dlatance/Inltlal  Momentum  Thickness,  x/0 

Fig.  25  Comparison  of  Prediction  with  Experiment,  NASA  Case  14, 
Two-Dimensional  Incompressible  Wake 

The  prediction  of  the  centerline  velocity  for  the  incompressible  planar  wake  is  shown 
in  Fig.  25,  plotted  in  similarity  coordinates.  From  similarity  considerations,  the  curve  of 
[1  -  (uc/ue)]'2  versus  x/0  should  be  a  straight  line,  and  Fig.  25  shows  that  the  prediction 
is  a  straight  line  in  these  coordinates.  However,  the  slope  of  this  curve  may  be  somewhat 
too  low,  indicating  too  slow  an  approach  to  the  true  asymptotic  condition.  This  effect 
was  also  seen  in  the  axisymmetric  jet,  and  it  is  interesting  to  note  that  the  prediction 
of  these  asymptotic  flows  made  by  this  model  are  quite  similar  to  the  predictions  of 
the  "k"  model  of  Launder,  et  al.  (Ref.  27)  as  presented  at  the  NASA  conference.  The 
"k"  model  of  these  authors  uses  an  eddy  viscosity-type  relationship  between  the  shear 
stress  and  the  kinetic  energy  and  an  algebraic  length  scale  formulation.  Thus  the  similar 
behavior  of  the  two  models  in  asymptotic  flow  indicates  that  the  length  scale  formulation 
may  be  at  fault;  i.e.,  it  may  be  necessary  to  use  a  more  sophisticated  formulation  including 
a  transport  equation  for  dissipation  (such  as  the  formulation  of  Hanjalic  and  Launder, 
Ref.  19)  in  order  to  properly  model  the  transition  to  asymptotic  flow. 

Figures  26  and  27  present  a  comparison  of  predicted  and  measured  velocity  and 
shear  stress  profiles  for  the  incompressible  two-dimentional  wake.  In  general,  the  velocity 
profile  prediction  is  good,  although  a  tendency  exists  to  underpredict  the  wake  width 
(which  is  consistent  with  the  underprediction  of  the  velocity  decay  rate).  Figure  27  shows 
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tlxat  the  prediction  of  the  peak  shear  stress  is  quite  reasonable  (and  it  should  be  remembered 
that  these  calculations  begin  with  estimated  rather  than  experimental  shear  stress  profiles), 
but  the  profile  width  is  again  underprcdicted. 


O 


Fig.  26  Comparison  of  Measured  and  Predicted  Velocity  Profiles, 

NASA  Case  14,  Two-Dimensional  Wake 

For  the  axisymmetric  incompressible  wake,  similarity  considerations  require  that  the 
curve  of  [l-(uc/ue)]"3/2  versus  x/D  be  a  straight  line,  and  Fig.  28  shows  that  such  a  curve 
is  predicted  by  this  model.  It  might,  however,  be  noted  that  asymptotic  conditions  for 
axisymmetric  wakes  have  not  been  experimentally  demonstrated  (Ref.  1),  and  thus  a 
judgment  as  to  the  proper  slope  of  the  decay  curve  cannot  be  made  from  the  evidence 
of  Fig.  28. 

Figure  29  shows  that  the  asymptotic  prediction  of  the  centerline  velocity  for  the 
compressible  two-dimensional  wake  again  shows  too  low  a  rate  of  increase  of  the  centerline 
velocity.  This  is  again  consistent  with  the  "k"  model  of  Ref.  27.  On  the  other  hand, 
the  prediction  of  the  axisymmetric  compressible  wake.  Fig.  30,  which,  as  discussed  above, 
may  not  be  an  asymptotic  flow,  is  considerably  better,  and  also  considerably  better  than 
the  "k"  model  of  Ref.  27. 
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Fig.  28  Comparison  of  Prediction  with  Experimental  Data,  NASA  Case  15, 
Axisymmetric  Incompressible  Wake 
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4.7  CONCLUSIONS  AND  COMPARISONS  WITH  OTHER  MODELS 

The  predictions  shown  above  have  demonstrated  that  the  present  model  is  capable 
of  accurate  predictions  of  a  variety  of  flows  of  engineering  interest.  Its  major  failures 
appear  to  be  in  the  prediction  of  the  asymptotic  circular  jet  and  the  asymptotic  circular 
wake.  In  this  respect  the  present  model  can  be  compared  to  the  similar  "k"  model  of 
Ref.  27,  but  the  present  model  provides  considerably  better  prediction  than  the  "k"  model 
of  such  flows  as  the  axisymmetric  wake  and  the  two-dimensional,  two-stream  jet.  This 
may  indicate  that  in  asymptotic  flows  the  algebraic  length  scale  formulation  is  in  error, 
since  in  such  weak  shear  flows  the  dissipation  term  in  the  turbulent  kinetic  energy  equation 
dominates.  On  the  other  hand,  the  improvement  of  the  present  model  over  the  "k"  model 
of  Ref.  27  in  strong  shear  flows  may  indicate  that  the  shear  stress  formulation  used  for 
this  model  is  more  accurate  than  the  eddy  viscosity  formulation  of  Ref.  27.  It  should 
also  be  noted  that  all  of  the  present  calculations  were  begun  from  estimated  initial 
conditions,  whereas  the  calculations  of  Ref.  27  were  begun  by  iterating  on  the  initial 
condition  until  the  data  at  the  first  downstream  station  were  matched. 

The  calculations  made  using  the  present  method  may  also  be  compared  to  the 
closely-related  integral  model  calculation  of  Ref.  12,  also  presented  at  the  NASA 
conference.  In  general,  the  integral  model  cannot  handle  the  developing  part  of  the  first 
regime,  in  which  boundary-layer-like  profiles  are  developing  into  profiles  characteristic  of 
a  free  shear  layer.  However,  in  flow  regions  where  the  assumption  of  cosine  profiles  used 
in  formulating  the  integral  model  is  satisfactory,  the  integral  model  predictions  are  as 
good  as  those  shown  here.  Further,  the  integral  model  prediction  of  the  asymptotic  jet 
(the  two-dimensional  wake  was  not  attempted  in  Ref.  12)  is  quite  good.  The  difference 
between  the  integral  calculation  and  the  present  model  in  this  flow  is  that  the  TKE  profiles 
assumed  in  the  integral  model  imply  a  much  greater  diffusion  rate  of  TKE  at  the  centerline 
than  is  predicted  in  the  finite-difference  model.  Thus,  although  a  comparison  of  the  present 
model  and  the  "k"  model  of  Ref.  27  put  the  length-scale  formation  for  an  asymptotic 
jet  under  suspicion,  a  comparison  of  the  present  model  with  the  integral  formulation 
indicates  that  the  proper  behavior  could  also  be  obtained  through  use  of  a  variable  TKE 
diffusion  coefficient. 

In  general,  however,  the  predictions  of  the  present  model  for  flows  of  engineering 
interest  are  quite  adequate;  they  are  certainly  better  than  those  obtainable  with  eddy 
viscosity  models  (Ref.  1).  The  present  model  is  also  adaptable  to  more  complex  cases, 
such  as  reacting  turbulent  flows,  using  both  conventional  and  unconventional  models  for 
the  appropriate  chemistry  (Refs.  28  and  29).  It  does  not  predict  asymptotic  flows 
adequately,  but  in  view  of  the  behavior  of  more  sophisticated  models  (Refs.  19  and  27) 
in  these  flows,  it  does  not  seem  worthwhile  to  attempt  to  improve  the  present  model 
for  asymptotic  predictions  which  are  of  limited  utility. 
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SECTipN  V 

THE  COMPUTER  PROGRAM 
5.1  MECHANICS  OF  THE  NUMERICAL  PROCEDURE 

The  sequence  of  operations  performed  by  the  computer  program  is  shown 
schematically  in  Fig.  31.  From  this  figure,  it  can  be  seen  that  the  operations  are  broken 
down  into  a  number  of  subroutines,  so  that  only  the  subroutines  involved  need  be  altered 
if  a  modification  is  desired.  In  general,  information  is  passed  between  the  subroutines 
through  the  use  of  labeled  common  blocks,  and  thus  some  care  in  making  modifications 
needs  to  be  exercised  to  avoid  altering  information  needed  elsewhere.  Table  IV  lists  the 
variables  stored  in  common  blocks  and  their  use  in  each  subroutine.  It  may-  be  noted 
from  this  table  that  several  variables  defined  in  the  original  program  (Ref.  13)  are  no 
longer  used,  but  the  appropriate  blocks  have  been  retained.  Little  effort  has  been  expended 
to  optimize  this  program  with  respect  to  run  time;  however,  a  typical  time  for  one  of 
the  computations  described  in  the  preceding  section  is  one  minute  on  an  IBM  370/155, 
and  the  program  occupies  86K  bytes  of  storage. 

The  sequence  of  computer  operations  for  a  typical  run  begins  with  those  operations 
that  are  only  performed  once  in  the  solution  of  a  particular  problem.  Referring  to  Fig. 
31,  the  appropriate  values  of  the  material  constants,  such  as  the  gas  constant  and  the 
specific  heat,  are  established  in  CONST.  The  heading  is  read,  and  all  other  operations 
for  which  no  specific  subroutine  name  is  shown  are  performed,  in  the  MAIN  program. 
This  heading,  which  may  consist  of  up  to  80  characters  (blanks  included),  is  printed  out 
at  the  top  of  each  page  of  output  for  identification. 

Subroutine  BEGIN  also  falls  into  the  category  of  subroutines  called  only  once  in 
a  given  problem.  In  this  subroutine  the  initial  profiles  of  velocity,  eddy  viscosity,  total 
enthalpy,  and  concentration  of  jet  species  are  read  in;  coefficients  are  initialized;  and  the 
lateral  grid  spacing. in  nondimensional  stream  function  (cu)  coordinates  is  computed,  using 
subroutine  DENSTY  to  obtain  the  densities  appropriate  to  the  initial  composition  and 
enthalpy. 

Length  scales;  the  mixing  region  width  YL,  defined  as  the  distance  between  the  points 
at  which  the  local  velocity  is  a  specified  function  of  each  edge  velocity  (see  Eq.  (17)); 
and  the  half-velocity  width  YMU,  defined  as  ihe  distance  from  the  centerline  to  the  point 
at  which  the  velocity  is  equal  to  the  average  of  the  free-stream  and  centerline  velocities, 
(see  Eq.  (19))  are  defined  in  subroutine  LENGTH.  A  third  length  scale,  based  on  an  assumed 
cosine  profile  shape,  will  also  be  computed  in  subroutine  SHEARS,  using  Eq.  (18).  These 
length  scales  are  necessary  for  the  computation  of  the  TKE  dissipation  rate,  which  is 
carried  out  in  subroutine  SOURCE. 
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Fig.  31  Computer  Program  Flow  Diagram 
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TABLE  IV 

VARIABLES  STORED  IN  COMMON  BLOCKS 
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TABLE  IV  (Continued) 


Block 

Variable 
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COEFF 

CONST 

| 

M 

1 

LENGTH 

MOMI2 

OUTPUT 

1 

SHEARS 

SLIP 

SOURCE 

VEFF 

HED 

HEAD 

I 

| 

1! 

■ 

I 

N 

I 

9 

NP1 

G 

U 

U 

9 

u 

U 

U 

NP2 

U 

G 

U 

U 

U 

U 

u 

U 

U 

NP3 

u 

G 

u 

U 

' 

u 

u 

U 

NEQ 

u 

I 

u 

IB 

U 

NPH 

u 

G 

u 

U 

KEX 

u 

I 

i.Rpj 

u 

U 

KIN 

u 

I 

liy! 

u 

u 

M 

13 

U 

KASE 

G 

■ 

KRAD 

u 

I 

1 

u 

u ! 

u 

U 

KPRAN 

I 

nl 

IDIN 

INDIC 

G 

|B 

IW 

KWAKE 

I 

■ 

1 

u 

ISL 

u 

I 

■ 

u 

KJU 

KMU 

■ 

G 

L 

AK 

ALMG 

isi 

LIME 

n 

□ 

IE 

M 

LI 

YL 

G 

U 

UMAX 

in 

G 

' 

UMIN 

u 

■ 

G 

u 

FR 

G 

LB 

U 

YIP 

IB 

G 

YEM 

>■ 

G 

OUT 

UP 

, 

G 

ifi 

LB 

u 

US 

G 

IB 

PHIUC 

G 

IB 

PSIU 

G 

HB. 

B 

u 

u 

XOD 

G 
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I 
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u 
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I 

RI 
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U 
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■ 
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■ 

■ 
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TABLE  IV  {Concluded) 


Legend:  G  -  Generated  in  subroutine 
M  -  Modified  in  subroutine 
U  -  Used  in  subroutine 
I  -  Input  data  to  subroutine 

NOTE:  Subroutine  SOLVE  does  not  use  common  block  storage 

After  the  length  scale  computation  the  turbulent  shear  stress  SHEAR  (I)  is  computed 
from  the  turbulent  kinetic  energy  F(1,I)  in  subroutine  SHEARS.  In  the  first  step  of  the 
computation,  an  additional  calculation  is  also  carried  out  in  SHEARS  in  order  first  to 
compute  the  value  of  SHEAR(I)  from  the  input  profiles  of  eddy  viscosity  and  velocity, 
and  then  to  compute  the  turbulent  kinetic  energy.  The  value  of  the  turbulent  Reynolds 
number,  RTM,  at  the  step,  which  is  used  in  computing  ASD2  and  the  TKE  dissipation, 
is  also  obtained  in  SHEARS. 

Subroutine  ENTRN  is  next  called  to  establish  the  viscous-region  mass  increase  over 
the  next  step.  As  explained  in  Section  III,  this  subroutine  is  fundamental  to  the 
computational  scheme,  yet  any  reasonable  method  of  obtaining  the  entrainment  rates  is 
acceptable. 

With  the  entrainment  for  the  next  step  computed,  the  step  size  is  obtained,  in  MAIN, 
by  demanding  that  the  step  size  be  such  that  the  mass  increase  in  the  viscous  region 
be  a  specified  fraction  (FRA)  of  the  viscous  region  mass  at  the  beginning  of  the  step. 
Because  this  length,  uncontrolled,  can  be  too  large  for  good  accuracy,  the  step  length 
is  externally  limited  in  MAIN  to  a  specified  fraction,  DXCN,  of  the  mixing  region  length 
scale. 


The  coefficients  necessary  for  a  solution  of  the  equations  at  X  =  XD  =  XU  +  DX 
are  then  obtained  in  COEFF  based  on  values  of  the  dependent  variables  at  X  =  XU. 
To  obtain  these  coefficients,  COEFF  calls  two  additional  subroutines,  SOURCE,  in  which 
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the  source  terms  in  the  transport  equations  are  computed  and  SLIP,  in  which  the 
coefficients  at  the  points  2  and  N  +  2  across  the  stream  are  evaluated.  Special  evaluations 
are  made  at  these  points  since  they  are  designed  to  produce  not  the  proper  values  of 
the  dependent  variables  at  2  and  N  +  2,  but  the  proper  gradients. 

This  point  marks  the  end  of  a  computational  step,  and  thus  profiles  of  the  dependent 
variables  and  some  of  the  coefficients  are  printed  out  through  subroutine  OUTPUT.  To 
start  the  calculations  at  the  next  step,  all  of  the  equations  involved  are  solved 
simultaneously  (except  for  the  momentum  equation  which  does  not  involve  other  variables) 
in  SOLVE.  Two  passes  through  SOLVE  are  made,  the  first  obtaining  the  downstream 
values  of  the  velocity,  and  the  second,  the  downstream  values  of  the  other  variables.  After 
this  computation,  the  calculated  profiles  are  inspected  for  possible  errors,  such  as  negative 
values  of  a  positive-definite  quantity,  and  any  errors  encountered  are  corrected. 

With  new  profiles  of  the  dependent  variables  known,  the  downstream  value  of  the 
viscous-layer  mass  flux  is  obtained.  On  the  first  step,  this  is  accomplished  by  summing 
the  contributions  of  the  entrained  mass  flux  and  the  mass  tlux  at  the  start  of  the  step. 
For  all  subsequent  steps  the  condition  that  the  integral  momentum  flux  increment  be 
a  constant  is  used,  which  involves  use  of  subroutine  MOM  12. 

The  values  of  the  physical  coordinate  Y(l)  which  correspond  to  the  nondimensional 
coordinates  OM(l)  and  the  computed  values  of  the  velocity,  total  enthalpy,  and 
concentration  are  now  established  in  subroutine  READY.  This  in  turn  involves  subroutines 
RAD,  in  which  the  radius  of  the  inner  edge  of  the  mixing  layer  is  obtained  and  DENSTY, 
in  which  the  density  RHOH)  corresponding  to  the  computed  values  of  velocity,  total 
enthalpy,  and  concentration  is  evaluated.  With  all  values  of  the  dependent  and  independent 
variables  known  at  the  new  X-station.  which  now  becomes  XU,  the  computation  is  repeated. 
The  sequence  is  continued  until  cither  the  maximum  number  of  steps  is  exceeded  or  until 
the  input  maximum  value  of  X(XL)  is  achieved. 

5.2  INPUT  DEFINITIONS 

All  initial  profile  data  are  read  in  in  either  the  main  program  or  subroutine  BEGIN. 
The  first  card  read,  from  MAIN,  defines  NCASE,  the  number  of  cases  to  be  computed 
on  this  run  of  the  program.  A  case  is  defined  as  a  complete  calculation  of  a  particular 
How  field.  NCASE  is  read  on  an  15  format.  The  second  card  supplies  the  heading,  with 
a  maximum  of  80  characters,  blanks  included.  Finally,  the  last  card  read  in  MAIN  supplies 
the  values  of  XOD1  and  URAT.  These  variables  are  used  in  the  case  when  a  calculation 
is  begun  at  some  x/D  downstream  (XODI)  with  a  nonunity  ratio  of  local  centerline  velocity 
to  jet  velocity.  The  variable  URAT  then  relates  the  computed  velocity  ratio  to  the  jet 
velocity.  These  are  read  on  a  2E12.5  format. 

The  remaining  data  needed  to  start  a  calculation  are  read  in  subroutine  BEGIN.  The 
first  of  these  is  a  set  of  10  integers  which  control  the  calculation,  read  on  a  1015  format. 
These  variables,  and  their  values,  are  as  follows- 
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TABLE  V 

INTEGER  CONTROL  INPUTS 


KRAD 

KIN 

KEX 

NEQ 


N 

^ISL 

k-' 

KPRAN 

KSST 

KSEV 

IWAKE 


=  0  plane  flow 
=  1  axisymmetric  flow 

=  1  inside  edge  of  flow  (I  *»  Da  wall  ^ 

=  2  inside  edge  of  flow  (Is  1)  e  free  stream 
=  3  Inside  edge  of  flow  (I  =  1)  an  axis  of  symmetry 

Similar  to  KIN,  but  applies  to  outer  edge  of  flow  (I  *  NP3) 
Number  of  equations  to  be  solved 

=  2  momentum  and  TKE> 

=  3  momentum.  TKE,  and  total  enthalpy 
c  4  momentum,  TKE,  total  enthalpy,  and  species 

Number  of  cross-stream  input  points  less  one 

*»  0  flo^v  is  a  2-D  shear  layer 
\  0  .any  other  fldw  field  type 

Not  used  m  presen:  formulation 

=  1  single -at  ream  flow,  US  ■  0 
X  1  no- effect 

=  0  UP  -  U(l) 

X  0  UP  =  UW 

*  1  flow  ls  a  wake 
X  l"any  other  flow 


After  these  variables  are  defined,  a  set  of  real  variable  controls  are  read.  The  definitions 
of  these  follow  in  Table  VI;  they  are  read  on  an  11F5.0  format. 


The  input  profij^s  of ^ Y( I ) ,  in  feet,  U(I),  in  ft/sec.  F(l,  I),  at  this  point  the  eddy 
viscosity,  in  f^/seCT^U,  !),  .the  total  enthalpy,  in  ft2 /sec2,  and  1  the  j*et  species 
concentration,  F(3,  I),  (dimensionless)  are  next  read  in  using  a  sequence  of  statements 
in  BEGIN.  They  are  read  sequentially,  using  a  7F10.5  format.  With  these  data  stored, 
the  calculation  begins. 


TABLE  VI 

REAL  VARIABLE  CONTROL  INPUTS 

XL.  Length  to  which  calculation  is  to  be  carried,  ft 

XPC6  Not  used 

ASD1  Value  of  ai,  normally  0.  3 

ASD2  Initial  value  of  a2.  normally  1.  69 

PREF(l)  TKE  Prandtl  number,  normally  0.  7 

PREF(2)  Mean  energy  trubulent  Prandtl  number,  0,7-0.85 

PREFO)  Turbulent  Schmidt  number,  normally  equal  to  PREF(2) 

DXC  Step  size  control,  generally  0.  5- 1.  0 

RO  Jet  nozzle  rsdius,  ft  —  'O  9 

HI  Inner  boundary- layer  thickness,  ft— 

CW  Characteristic  velocity,  normally  nozzle  exit  velocity 

5.3  VARIABLE  NAMES  AND  LISTING 


A  complete  list  of  the  FORTRAN  variable  names  used  in  the  program  and  their 
meanings  follows  as  Table  VII.  Appendix  I  provides  a  complete  FORTRAN  listing  of  the 
program  used  to  compute  all  of  the  test  cases  of  Section  IV. 
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TABLE  VII 

FORTRAN  VARIABLE  NAMES 


FORTRAN 

Name 

Mathematical 

Equivalent 

Description 

Defining 

Subroutine 

A(J,1) 

A 

Coefficient  m  difference  equation 

COEFF 

*ADUDY(D 

|du/8y! 

Absolute  value  of  velocity  gradient 

SHEARS 

•ADUDYM 

|3u/3y|  M 

Maximum  value  of  I3u/3y  at  x-8tation 

SHEARS 

♦AM 

u  -  u^Jdu 

Momentum  wcrement/maBS  flux  at  x-station 

MOMI2 

♦AMCHK 

<V,E  ■  ♦IV J,u  '  ue)du 

Momentum  increment  at  x- station 

MAIN 

♦AMD 

--- 

Value  of  AM  at  XD 

MAIN 

AME 

“E 

Mass  flux/ unit  area  entrained  in  AX  at  external 
boundary 

ENTRN 

AMI 

*1* 

As  above,  but  internal  boundary 

ENTRN 

♦AMUP 

... 

Value  of  AM  at  XU 

MAIN 

♦ASD1 

al 

Parameter  in  t  -  k  relation  ,  Eq.  (9) 

SHEARS 

+ASD2 

“2 

Parameter  m  TKE  dissipation,  Eq.  (14) 

SHEARS 

AU(I) 

Au 

Coefficient  in  velocity  difference  equation 

COEFF 

w.n 

B 

Coefficient  in  difference  equation 

COEFF 

BETA 

... 

Nut  used 

BEGIN 

BU(D 

Bu 

Coefficient  in  velocity  difference  equation 

COEFF 

*c 

1 

C1 

Density  correction  to  TKE  dissipation 
(Eqs.  (23  and  24)) 

SHEARS 

C(J,I) 

C 

Coefficient  in  difference  equation 

COEFF 

CPF 

cPf 

Frozen  value  of  constant  pressure  specific 
heat 

DENSTY 

CPI 

... 

Jet  specific  heat 

CONST 

CP2 

... 

Free-stream  specific  heat 

CONST 

GSALFA 

Cos  a 

Cosine  of  angle  between  flow  and  axis  of 
symmetry 

RAD 

CU 

.  Cu 

Coefficient  m  velocity  difference  equation 

COEFF 

DEN 

P  ret 

Reference  value  of  density,  not  used 

... 

DPDX 

dP/dx 

Pressure  gradient,  here  ■  0 

MAIN 

DUDOM 

Bu/Bw 

Velocity  gradient  in  i«;  coordinates 

MAIN 

*DUDY  ' 

du/8y 

Velocity  gradient,  ft/ sec /ft 

SHEARS 

DX 

AX 

Axial  step,  ft 

MAIN 

*DXC 

... 

Max  step  length  factor 

MAIN 

*DXCN 

... 

DXC/Z 

MAIN 

*DX1 

... 

Value  of  DX  before  limiting 

MAIN 

*EPSI 

eI 

Constant  eddy  viscosity  near  inside  edge 

SHEARS 

•EPSO 

CE 

Constant  eddy  viscosity  near  outside  edge 

SHEARS 

EMU 

PT 

Turbulent  dynamic  viscosity 

VEFF 

F(1,I) 

K 

Turbulent  kinetic  energy,  k,  (ft/sec)2 

SOLVE 

F(2, 1) 

H 

Total  enthalpy  (ft/sec)2 

SOLVE 

F(3, 1) 

C 

Jet  species  concentration,  dimensionless 

SOLVE 

FRA 

... 

Constant  m  DX  expression 

MAIN 
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TABLE  VII  (Continued) 


FORTRAN 

Name 

Mathematical 

Equivalent 

Description 

Defining 

Subroutine 

GCl 

Hi 

Gas  constant  for  jet 

CONST 

GC2 

R2 

Gas  constant  for  free  stream 

CONST 

♦HEAD 

... 

Heading  {20 A4)  format 

MAIN 

*IE 

— 

Outer  edge  point  at  which  entrainment  is 
calculated 

ENTRN 

♦IECSE 

... 

Print  control 

MAIN 

*11 

... 

As  IE, 

for  inside  edge 

ENTRN 

♦IK  MAX 

... 

Point  In  profile  where  t  -  Tmax 

SHEARS 

♦1MAX 

. 

Point  in  profile  where  DUDY  =  DUDYmax 

SHEARS 

♦IND1C 

... 

Determines  case  number 

MAIN 

INTG 

— 

Axial  step  number 

MAIN 

♦IPHNT 

... 

Print  control  step  counter 

MAIN 

♦1SL 

... 

Shear 

layer  indicator 

BEGIN 

♦IWAKE 

... 

Wake  indicator 

KASE 

... 

Free- 

or  wall -flow  indicator 

KEX 

... 

External  boundary  designator 

KIN 

... 

Internal  boundary  designator 

♦KPRAN 

... 

Not  used 

KRAD 

... 

2-D  or  axisymmetric  indicator 

*KSEV 

... 

Characteristic  velocity  control 

,  *KSST 

... 

Single 

-stream  Indicator 

N 

... 

Number  of  increments  across  stream 

1 

♦NCASE 

... 

Number  of  flows  to  be  computed 

MAIN 

NEQ 

... 

Number  of  transport  equations  to  solve 

BEGIN 

NPH 

... 

NEQ-1 

BEGIN 

NP1 

... 

N-r  1 

NP2 

... 

N  +  2 

BEGIN 

NP3 

... 

N  +  3 

OM(l) 

<*  -  -  #i> 

Independent  cross -stream  coordinate  {non- 
dimensional  stream  functicn) 

BEGIN 

OMD 

Ui+I  '  ui-l 

Increment  in  u 

COEFF 

•  P(l) 

... 

Not  used 

... 

PEI 

^E  "  *1 

Mass  flux  in  mixing  region 

MAIN 

♦PEI1 

--- 

PEI  at  XU 

MAIN 

♦PE  12 

... 

PEI  at  XD 

MAIN 

♦FHIUC 

(u  -  u  )/(u  -  u  ) 
c  e  j  e 

Nondimens ional  centerline  velocity  ratio 

MAIN 

PRO) 

... 

Not  used 

... 

PREF(l) 

Prk 

TKE  Prandtl  number 

BEGIN 

PREF(2) 

PrT 

Turbulent  Prandtl  number  (mean  energy) 

BEGIN 

PREF<3) 

ScT 

Turbulent  Schmidt  number 

BEGIN 

♦  PS  1IE 

tnB!-V/ta  c'V 

Nondimens  ional  velocity  at  outer  entrainment 
point 

ENT.HN 
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TABLE  VII  (Continued) 


FORTRAN 

Name 

Mathematical 

Equivalent 

Description 

Defining 

Subroutine 

•psin 

(un  -  ue>/(uc  -  ue> 

As  PSILE,  at  inner  entrainment  point 

ENTRN 

•PS1U 

(u  -  u  )/(u  -  u  ) 

e  c  e 

Nondimensional  velocity 

MAIN 

•PS1UD 

(a  -  u  >/(u  -  u  j 
e  c  e 

Desired  nondimensional  velocity  at  given  en¬ 
trainment  point 

ENTRN 

RID 

r 

Radial  coordinate  »  RI  Y(I)*CSALFA 

READY 

•RAAUHfD 

r2  p  u 

... 

SHEARS 

RBAR 

R 

Mixture  gas  constant 

DENSTY 

*RD 

rD 

Inner  edge  radius  at  XD 

MAIN 

RHOU) 

P 

Density  Gbm/ft3) 

DENSTY 

RI 

rI 

Radius  to  inner  edge  of  mixing  region,  ft 

RAD  <7- 

HO 

r 

o 

Nozzle  radius,  ft 

BEGIN 

•RR1 

... 

Turbulent  Reynolds  number 

SHEARS 

»RU 

r 

u 

Inner  edge  radius  at  XU 

MAIN 

SC(I) 

c 

Coefficient  in  diffusion  term  of  general  equa¬ 
tion 

ENTRN 

•SCSHU) 

... 

Modified  coefficient 

ENTRN 

•SHEAR  (I) 

T 

Turbulent  shear  stress 

SHEARS 

•SIGMA 

a 

2-D  shear  layer  growth  rate 

SHEARS 

TEMPID 

t 

Local  static  temperature 

DENSTY 

•TOJO 

T  ./T 
oj  oe 

Centerline/free-stream  total  temperature 
ratio 

OUTPUT 

U(I) 

u 

Velocity,  ft/ sec 

SOLVE 

•UDES 

a 

Value  of  u  desired  at  given  lateral  point  at 
downstream  stations 

ENTRN 

•UGD 

... 

Free -stream  velocity  at  XD 

MAIN 

•UGU 

... 

Free -stream  velocity  at  XU 

MAIN 

•UMAX 

... 

Maximum  velocity  in  a  given  profile 

LENGTH 

•UMIN 

... 

Minimum  velocity  in  a  given  profile 

LENGTH 

•UP 

U. 

J 

Centerline  velocity  at  X  =  0 

BEGIN 

•URAT 

Coefficient  modifying  PHIUC 

MAIN 

•US 

u 

e 

Secondary  stream  velocity  at  X  =  0 

BEGIN 

•UTOL 

... 

Minimum  velocity  difference  at  a  step 

LENGTH 

*uw 

... 

Reference  velocity 

BEGIN 

X 

X 

Axial  position,  ft 

MAIN 

XD 

*D 

Axial  position  at  downstream  station 

MAIN 

XL 

L 

Value  of  x  at  vehicle  to  stop  calculation 

BEGIN 

•XOD 

x/D 

Nondimensional  axial  coordinate 

MAIN 

•XODI 

x/Di 

Initial  value  of  x/D 

MAIN 

XP 

... 

X  at  last  step 

MAIN 

•XPC6 

... 

Not  used 

i 

_ 1Z _ 1 
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TABLE  VII  (Concluded) 


FORTRAN 

Name 

Mathematical 

Equivalent 

Description 

Defining 

Subroutine 

*XSTA 

... 

X  at  which  eddy  viscosity  profile  modifications 
are  used 

SHEARS 

XU 

X 

u 

Axial  position  at  upstream  station 

MAIN 

*XXU 

... 

XU,  in. 

MAIN 

Y(I) 

y 

Physical  lateral  coordinate,  ft 

READY 

*YD 

yD 

Value  of  YI  at  downstream  step 

MAIN 

*YI 

yi 

Distance  to  inner  edge  of  2-D  mixing  region 

RAD 

*YL 

i 

2-D  region  width  scale,  Eq.  (17) 

LENGTH 

*YLK 

Dissipation  length  scale,  Eq.  (17  through  19) 

SHEAR 

*YMU 

y  1/2 

Value  of  y  at  which  u  =  1/2  (u  +  u  ) 
c  e 

LENGTH 

*YY(D 

... 

As  Y(I),  in. 

MAIN 

*YU 

... 

Value  of  YI  at  upstream  step 

MAIN 

indicates  variable  not  used  in  original  program  (Ref.  13) 
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86  G  =  IU<  NP2 I »U( NP1 1-8. *U (NP3 I  1/ (5. * I U (NP 2 1 *U INP1 )) *8.*U ( NP 3 ) ) 

0089 

GF  =  I 1 ,-PREF( JI)/(1.*PREF( J)  1 

0090 

T7F=  r(T*GF  1 71 1 .  *Tr*GF1 

0091 

K J,NP2)=F( J.NPl l«GF+(l.-GF)«F( J.NP3I 

0092 

GO  TO  88 

009  3 

87  F(  J,NP2I«(4.*F1J,NP31-F  (J.NPDI/3. 

0094 

88  CONTINUE 

0075 

45  CONTINUE 

0W5 

■(TSLT  UCNSTY  ' 

C 

CALCULATION  OF  RAOI  ( 

0097 

CALL  RADIX U,R(1),CSAL FA) 

0098 

IFICSALFA.EO. 0. .OR.KRAC.EG.O)  GO  TO  27 

0099 

00  28  1=2, NP3 

0100 

28  R I I 1 *R( l !♦ VI I )*CSALF A 

trial 

GO  TO  29 

0102 

27  00  30  1=2, NP3 

0103 

30  P II )=R(1) 

0104 

29  CONTINUE 

.  C 

CALCULATION  OF  OMEGA  VALUES 

0105 

OMI 11=0. 

0106 

OMI  21  =  0. 

0107 

DO  49'  I  =3  ,NP2 

0T3B 

49  OMI  U=UH1  [-1 J+.  5*(RH0II  J+U(  l  l*RI  I)  *RHOI  I- 1)*U'(  I-n*R(  1-1)  1* 

l(Y(Il-Y((-in 

0109 

PF I =0M(NP2) — 

0110 

00  59  ! =3 , NP1 

0111 

59  OMI 1 1 =0M( 1 ) /PEI 

0112 

OM ( NP2 I =1.0 

0113 

0HINP3) =1. 

0114 

RETURN 

0115 

END 
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<5001 
°002_ 
"0d03  ‘ 


SUBROUTINE  COEFF  . 

IMPLICIT  REAL*8!A-H*0-Z) 

COMMON  /GEN/PEI (AMI  ,AME , DPDX, PREF I3t » PRI3J, P 1 3) ,OEN, AMU.XU.XD, XP • 
lXL«OXf CSALFA t XPC6 ►  I NTG 


l/V/UI 43 ) ,F(3,43),RI43),RH0I43)»0HI43),YI43) 

1 /C/SC  1 43) . AUl 431 ,BU(43) ,CUI43 ). A( 3,43 ) ,8 ( 3,43) ,C! 3,43) 
COMMON  /SHEAR/  SHEAR  143) .SCSHI43) 

'  '  C0MM0N/DUD/DUD0MI43I ,  DUDYI43),  ADU0YI43),  AOUDYM 

DIMENSION  G 1 1 43 ) ,G2(43) ,G3(43),D( 3,43 )«S 1143) ,S2I43) , S31 43) 


THE  CONVECTION  TERM 
~  SA=R{ 1) *AMI/PEI 

SB*I RINP3 )• AME-RI1 ) *AMI )/PEI 
'  OX=XD-XU 

00  71  I =3f NP1 


P2=. 25/OX 
Pi-92 /OMO 

Pl-U)MtI  +  ll-0M(III*P3 
P3  =  I0M(I) -OH 1 1-1 ) ) *P3 
P2=3. *P2 


U=SA/UI”U 

R2»-SB*.25 

R3-R2/CMD 

Al*-l0MIIU>O.*CN(I) >*R3 
R3=lOM|I-l)F3.*CM|I  ))  *R3 
GUI  )=P1*U»R1 


v»^t  i  i 

G3II )=P3-Q+R3 

Cut I)=— Pl*UII*l)-P2*UII)-P3*UtI-l) 
THE  DIFFUSION  TERM 
AUID-2./0M0 

Bum=SCII-ll*AU(  I)/ICMII)-CMU-1)) 


MHiin.i«f  inii«i>fa:nmDu:n 


IFINEC.EQ.l)  GO  TO  33 
00  34  J-l.NPH 

Cl J.I )«-PI*F! J,l*l)-P2*FI J,I)-P3«FIJ,I-11 
CALL  SOUKCEICS,OIJ,II,J,I) 

Cl J.I)*-CIJ,I )*CS— F  f J (I)TOlj. I) 


11/  rntr  1  J  I 

BIJ,I)»BUII)/PREFI J) 

34  CONTINUE 

SOURCE  TERM  FOR  VELOCITY  EQUATION 
33  PHI  =  0.0 

Sill)  =  IOPOX  ♦  PHI )*0X 


L»JMTT>11 


S3II)«P3*S1I I ) / I RHO II -1 ) *UI I-l) > 

sitt  >cPi*sim/(RHoim)Mji!»i) ) 

CUII l*-CUII>-2.*ISltI>»S2II>*S3tI>> 
Sit I)*SII  II/UIIU) 

S21 1 1 =  S2 II) /UI I ) 


71  CONTINUE 

COEFFICIENTS  IN  THE  FINAL  FORM 
DO  91  I =3»NP1 

RL=l./(G2t I )*AUI l)*BU( I J-S2I I )) 
AUI I ) *( AUt I ) +S1 I I )-Gl I I ) ) *RL 
BUII)’1BU1I)  +  S31I )-G3 1 1 ))*RL 
91  CUI I )=CUI I )*RL 


I 

DO  92  J=1,NPH 
00  92  I-3.NP1 

RL-1./IG2ID  +  AI  J.Il+BIJ.I  >-0(J«I)  ) 
AlJtI)=IAIJ,I)-GllI))*RL 
BIJ,I)=IBI JtI)-G3II))*RL 


LI  Jtl l=LI 

76  CALL  SLIP 
RETURN 
ENO 


60 


AEDC-TR-73-1 77 


SLBROUTINE  CONST 
IMPLICIT  REAL*8(A-H,0-Z> 

COMMON  /GEN/PEI .AMI ,APE.DPDX7PftEF(3TTPR( 3T7P T3I ,'DEN, AMU, XU. XO,XP, 
1XL.DX,CSALFA,XPC6.INTG 


.inrliNrciNrJf NtW  t  nrn*  KcA  t 

1/Ll /VL, UMAX, UMIN,FR,YIP,Y EM 
COMMON  /ASD/AS01.AS02  # 

COMMON/C PS/CP1(CP2  «GC1»GC2  — - 

CCMMON/IOIN/INDIC  ® 

FR-.IO 

„*&■»■*>>■***  w*4* 

-  PA  ( 2 1  *  •  7  .  I ' 

AMU  =  0.000012  "fbr<f 

CP1*=3.42  - _ 

C P2=0. 24  —  ~  Q.'ZA Z. 

GC  1—766*6  - —  7  1  ^ 


GC2*53^35  -  N 

C  “RETURN  n?  T)$»  D*) 

END  L  ICPHUT  (  '  J  ' 

Pz^i£-M/lT 


6J 
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0001 

0002 

0003 


0004 

0005 

0006 

0007 


"ffdoTT 

0009 

0010 

0011 

0012 

0013 


0014 

0015 

0016 

0017 

0018 

0019 


0020 

0021 

0022 

0023 

0024 

0025 


0026 

0027 

0028 

0029 
00  30 
0031 
0032 
0033 
0034 
0035 
0036 


0077 
0038 
0039 
0040 
0041 
0042 
004  3 


SUBROUTINE  DENSTV 
IMPLICIT  RE AL*8 I A-H ,0-2 ) 

COMMON  /GEN/ PE I  , AMI  ,AMb , DPOX , PREF (31, PK 13), PI  31 , DEN, AMU , XU,  XD,XP , 
1XL,DX,CSALFA,XPC6, I  NT  G 

T/ V / U (  4 IT ,“F  I  3  , 4 3T r 4 3  r,  Rfl CT 4 3 )  ,  C M { 4 3 )  ,7143) 

1/I/N,NP1,NP2 ,NP3  ,NEU, NPH, KEX,KI N, KAS E ,KR AD, KPRAN 
CCMMON/A JXP/TEMPI43 ) 

COMMON/TEM/TOJO 
CCMM0N/CPS/CP1 ,CP2 , GC1 , GC2 
CCMMON/IW/IHAKE _ 

COMMCN/OENC/C 
DIMENSION  RBARI43I 
PINF=14. 7*144.  -v. 

00  45  1=1 , NP3 
_jLE_JNPH.LT.  3)  GO  TO  46 


~  ^  cff=c pr*i  i a n  i ♦ cps * 1 1  .-in ,  i  n — ■ 

CPF=CPF *2 5000.0 

R8AR I  II =GC1*F (3,1 )*GC2*(1.-FI3»  1) ) 

GO  TO  44 

46  CPF=. 24*25000. 

FI  3, 1 1=1. 

- p*AR(  [5=53.35 . . — 

IFINPH.LT. 2)  F(2,tl-CPF*520. 

44  TEMPI  I  )  =  (  F(2 , 1 ) -.5*1)1 1 )  *UI  I ) )  /CPF 
PHOII )=PINF/ITEMP11 )*«BARII ) ) 

IF(I.EO.l)  T1-FI2, 1)/CPF 

_ LH  1.EQ.NP3I  T2*F 12, NP3 l/CPF 

~45~X(TfrTTKUE - - ■  — 

T0J0=Tl/T2 

IFI I NTG.NE.O)  GO  TO  50 
CHECK  THAT  GAS  CONSTANTS  CORRESPOND  PROPERLY. 
I  J=1 
IA*=NP3 

- IFTUI II.GE.UINP3I IGO_Tir  4T 

I J=NP3 

I  A=1 

47  CONTINUE 

IFIRHOII J).LE.RHO(IA) ) 

TPAT=T2/T1 

- - 1FIIJ.E0.NP3)'  TRAT=T1/T2 - 

C  =0.95+. 05*K3AR (IA)*TRAT/RBARIIJ) 

GO  TO  50 

48  C  =0. 984+. 016*RH0(1A)/RHC( IJ1 
50  CONTINUE 

RETURN 

- nro - 


xruul 

.pCZ.eq.  1/  Tl5  TEHP&). 

Vf . 

^O07l^c?v 

?HoaV‘  Yi  * ?  rz)  J 

.  Lt,  O 

a(*  ctt* 


2=L0U  SPEED  SIDE 

r  c4>L 

“r 


4* 


GO  TO  48 


rU 
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OOJI 

0002 

0003 

0004 


0005 

TJUJV 

0007 

0004 

0009 

0010 

0011 

7TOTT 

0013 

0014 

0015 

001b 

0017 

MI5T 

0019 

0020 

0021 

0022 

0023 

THT2V 

0025 


0026 
7T07T 
0028 
0029 
00  10 
0031 
0032 


0033 

0034 

0035 

U03b 

”0037 - 

0033 

0039 

0040 

0041 

0042 

~0O45~  * 
0044 
0045 
0046 


SUBROUTINE  ENTRN 
IMPLICIT  REAL*8(A-H,C-Z> 

RE AL*4  HEAD 

COMMON  /GfcN/PEI  .AMl'.AME.DPDX.  PREF  (3 » , PR t 3) ,P < 3) .OEM, AMU, XU, XD, XP , 
1XL,DX,CSALFA,XPC6,INTG 

1/V/U14  3) ,F(3,43) ,R [43 1 , RHCI43 1 ,QH (43 1 ,Y(43) 

1 / 1 /N , N P I , NP2 ,  NP3 , NEQ. NPH.KEX , KIN.KAS E.KRAO,  KPRAN 
1/C/SC ( 43) tDUMM(5I61 
1/L1/YL  .UMAX.UM1  N.FK.Y'I  P.Y EM 
COMMON  /SHEAR/  SHE AR 143) .SCSHI43) 

CChMON  /RUH/  RAAUHI43) 

COMMON/LIME/ 1 1, IE 
COMMDN/PRNT/1 PRNT 
COMMON  /ASO/  AS0I.ASD2- 
COMMON/HEO/HEADI20) 

_ DO  99  *=2»NP2  _ 

SC  1 1 )=RAAUH( II*EMU/(PEI*PE1) 

99  CONTINUE 
SCI  11-0.040 
SCI NP3) *0.D+0 
IF ( IPRNT . EQ.01  GO  TO  10 
- WWTTET6T1  OOrrFTAD - 

1001  FORMA  Tin  DATA  FROM  ENTRN  ’  ,  20A4 1 
hRITF (6,1002) XI N.XEX, KPRAN, NEC, XU 

1002  FORMAT! >  Kl N=', )2,'KEX=» , 12, »KPKAN=» , 12, 'NE0=*, 12  ,  * XU«* , F8. 5 1 
10  CONTI  NIJF 

GO  TO  1 11  ,12  ,31  )  •  KI  N 

- ITGTT7ITV0 - - 

C 

12  IFIXU.EO.O.)  GO  TO  30 


UFR=I.01 

— rFruTii/u(NP3T.Lf.o;i)  ufr-i-.ict 
00  25  I*4,N 
I  1  =  ) 

1FIUA0SIUM l-UI 141J ) .GT..0001*UMAX)  GO  TO  26 

25  CONTINUE 

26  CONTINUE 


IFIUI  in.LT.UIlllGO  TC  251 
1FIUI Ill.GT.UIlllGO  TC  252 
IFIUI Tl4l).Lr.UllI)>GC  TO  251 
>F(U(  II  +  1).GT.U(IIMGC  TO  252 
C  THE  VALUE  OF  DOES  MAY  QE. VARIED  TC  SUIT  PROBLEM 

'  251  UOE S=0.  99"  *U(  1 )  . 

PSI UD=UDES 

IFIUmi.GT.UDES)  UDE  S=U  t  III 
GO  TO  253 

252  UOES=UFR*U! 1) 

1 F I UOES.EQ.O. )  U0ES-0.9*U(II ) 

- (>5TU0=U0ES 

IFIUI 1 1 1 . LT. UOt S )  UDES-UII1) 

253  CONTINUE 

226  GI  =  I2.*PE  I/IUIII+ll-UIII-l)))*!  I  (SCI  II  >*{U(  t  I*U-u(  I  m/(ON(  (  J  +  l 
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0047 

0048 

“00 49 - 

0050 
“0051 
0052 
17053 — 
0054 
"0055 — 
0056 

0057 

0058 

0059 — 

0060 — 

0061 — 


0062 

0063 

0064 

0065 

UU66 

0057' 

0068 

"0059 

0070 

0071“ 

0072 

0073 

0074 

0075 

0076" 

0077 

0078 

0079 

0030 

0081 

0032 

0083 

0084 

0085 


0086 

0087 


i)-0ft(rmi-(scm-i)*<uuii-uu  i-u  i  /  to*  t  in-on  m-im  i-ioMtii+i) 

1-0M(I I  —  1 1 t*  HIDES  -Util  1 )/IDX*2. 1 > 

Gt-DABSIGl) 

PSIUO=tPSIUO-Ut NP3I)/ (U (L ) — U ( NP3 I ) 
t>SI  1 1  =  1 U 1 1  I  I  -UI  NP3)  )  /  I U  ( 1  »-Ul5f*3i  ) 

GRAT-I 1.-PSIIII/I1.-PSIUDI 
IF(GRAT.GT.1..0R.GRAT.LE.0.)  GRAT=L. 

GI>GI*GRAT 

- IFTGl.GT.50.)  Gl-50. 

[Fill. GT. 51  Gl  =  0. 1*01 

- m'PRNTi'EO. 0)  GO'TC  40 - 

400  KRI TE 16,1000111  ,UI 1 1- 1 1 ,U 1 1  It ,U ( 1 1  +  1 1 ,0* <  II-l), OM (II 1 ,0M I 1 1  1 }  •  SC  < 

— rn=T)iscim.scMi*n,Gi 

1000  FORMAT!*  IN  ORDER  II,  U (31 ,CM(3 » , SC  I  31 , G I •/ 1 3, 10D11.4 J 

- GtnCTAOi  '  ‘ 

C 

- 30  "AWT^BAH5tT5HESRTZT*5HEinU5T;=2VS5H  E  A JET TT77 - 

I  _ I  U(2  J+U!3!-2.*U(  1)1 1 

GO  ro  40" 

c 

'  31  AHr-O. - 

c 

r 

40  GO  TO  (41 ,42 • 75  I , KEX 

41  IF  TKPRAN.EQ. l.QR. NEQ. EQ. 1 • OR.XU. EQ.O. )G0  TO  43 
IFIKIN.EQ.2) AM1=(GI-0PII I )*R( AP3 I  PANE )/  I l.O-OHI I I)«R(11) 

43  GO  TO  80' 

C 

- 42  TFTXU .TDTTJTT'TjCT  TO  70 - 

C 

00  56  1-3, N 

J-NP3-I 

IE-J 

I F I  DABS  I UI J )-U( •)♦  1 )  I  .GT .. 0001 -UMAX)  GO  TO  57 

- 56  CONTINUE" - 

57  CONTINUE 
C 

[Ft  U(  IE) •  LT .UINP3) ) GO  TO  551 
IF(U[IE).GT.U(NP3))G0  TO  552 
1  Ft  UI  IE-11  .LT.UdEDGC  TO  551 
IFlUt  IE-U.GT.UIIE)  IG0  TO  552 
C  THE  VALUE  UF  UOES  NAY  BE  VARIEO  TO  SUIT  PROBLEM 
551  UDES=0.99*U|NP31 
PSIUO-UDES 

I F (Ut I  El .GT.UDESI  UOES-UIIEI 
GO  TO  553 

— 552  UPC 5=1.0 1  *UINT>3I - 

PSIUD=UOES 

IFIUIIEI .LT.UOFSI  UOES-UtlE) 

553  CONTINUE 

556  GE- (2. *PEI/ lull E+l l-U ( 1 E- II I )*( ( f SCt IE  )*IUl IE*1 )-UI IE  11/ 

I I  OMf  !£♦!  )-  OMIlEm-(SC(lE-l)«tUUE)-UIIE-l))/ 

- ITTtHt  IE  I  -  OHIIE-im  1-tCHI  IE  +  1 J-OHI  IE-HJ*IUUES - -UIIEJT 

1/(0X*2.) I 
GE— l.*DA8S(GE) 

P  SI UO-t  PSI UO-UI NP3 1 )/ IU ( 1 1— Ut  KP3 1 1 
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0088 

0089 

0090 

oon 


0092 

0093 

0094 

0095 

0096 


P  SI  I  E  sl  U! IE)-UINP3> )/ IUI1)-UINP3> ) 

GRA T=PSI I E/PSI UD 

IFIGRAT.GT. 1..0R.GRAT.LE.0. I  GRAT=1. 

GE*GF  *GRAT 

lRdt.LT. INP3-5 ) )  Gt=0.1*GE 
IFIGE.LT.-100.)  GE»-100. 

I F { IPRNT.EQ.O)  GO  TO  60 

500  WRI TE(6tl003)IE,UIIE-l),U(IE)tU(IE+l) ,0M (I  E- 1  ) . OH ( I E I , UM ( IE *1 )  t  SC ( 
HE— li  > SC  I  IE)  .SCI  IE*1)  ,GE 

1003  FORMAT!*  (N  ORDER  IEt  U 13 1 , CMI3 > ,SC 1 3 1. GE* / 1 3, 1001 1 .4 1 


0097 

0098 

0099 

0100 


0101" 

0102 

0103 

0104 

0105 

0106 

orrr 


60  IF!KIN.EQ.1IAME  =  (GEHGM!IE|-1.)*RI1>*AMI  I/IRINP3)*  OMHEII 
1FIKIN.E0.2IG0  TO  65 
IFIKIN.E0.3I AME-GE/ !R (KP3 ) *GM ( 1 E ) ) 

'GO  TO  80 

C  80TH  BOUNDARIES  ARE  FREE 

- 65  AMI'ifGT *0H!  I  F1-GE*0HI  mi /TTCrFnn-UHin  iTTRTlTl - 

AME  =  ( { 1.-QH1 I  I  I )*GE-( l.-OMI IEI l*GI 1/ ( (0M( IEI-0M1 1  I ) |*R(NP3I  > 
lFdRRNT.ECi.OI  GO  TO  80 
WRITE (6. 1004) AMI ,AME,PEI 

1004  FORMATC'" AT  65  AMI  =  ■ ,D1 1 . 4, • APE*’ , Dl 1 .4, • PE  I* * t Dll. 4) 

GO  TO  80  / 

- Tff  AWE--DAB5n  Sh EARTKP21  *SHEJRT MWT=2'7*5HF«R  IDP3TT7 - ~ - 

1  !U1NP2>*UINP1)-2.*UINP3I> I 

IF ( I PRNT. EO. 01  GO  TO  80  ' 

702  WRI  TE 16.100 5) AMI .AME  « PE  I 

1005  FORMAT! •  AT  70  AMI** ,D1 1.4, • APE- • , Dl 1.4. • PEI* • , Dll. 4) 

GO  TO  80 


0108 

0109 

0110 

0111 


UTT2 - 

75 

AHE=0. 

0113 

IF! XU.EQ.O. ) 

GO 

TO  80 

0114 

IF! KI N.EQ.2) 

AMI 

■Cl/l!l. 

0115 

80 

RETURN 

xrn'6 

END  ~ 
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0001  SUBROUTINE  LENGTH 

0002  IMPLICIT  REAL*8I A-H.O-Z ) 

C003  COMMON  /GEN/ PE  1 , AMI ,AME, DPOX, PREE (3) , PR13J ,PI 3» , DEN, AMU, XU.XDt XP, 

_  1XL,DX,CSALFA,XPC6,INTG 


0031  CO  21  K*3,NP1 

0332  UMUZ=DABS (Ul KI-UHI 

0033  IFIUMUZ.LT.UMUZZ]  UMUZZ*=UMUZ 

0034  21  CONTINUE 

0035  UMU=UM*UMUZZ 


20  DO  22  K*3,NP1 

IFlUlKJ.NE.UHUl  GO  TO  22 
GO  TO  24 
22  CONTINUE 

UPU=UM-UHUZ2 


0037 

0038 

0039 

0040 

0041 


0043  IFtKRAZY.EU.2l  GO  TO  20 

0044  24  KKU=K 

0045  IF(UIKKU).EO.UM)  Y  MU=Y  ( KKU) 

0046  1FIUI 1I.GT.UINP3) J  GO  TO  25 

0047  I  FI Ul  KKUI.GT.UMI  GO  TO  28 

TO41 - IFIUIKKUI . LT. UH)  GU  TC'  27 - 

0049  25  IFIUIKKUI.GT.UMI  GO  TC  27 

0050  IFIUIKKUI.LT. UM)  GO  TC  28 

0051  27  YMU^YIKKU)*! UM-UI KKUI I *1 Y IKKU^l l-YIKKUJ ) / (UIXKU+I I-UIKKU) » 
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0052  GO  TO  61 

0053  28  YMU=Y  (  KKUI  +  I UM-LIl  KKU)  I  *{Y  l  KXU-l )  —V  ( KKU)  1/ lUtKKU-l»-UlKKl)>  ) 

0054  61  CONTINUE 

0055  DIF=DABSI UMAX-UMI N) *F R 


r 

0056 

IFIKIN.NE.2)  GU  TO  43 

0057 

47  J-2 

0058 

48  J=J*1 

0059 

UJ1*UI JI-UIl) 

0060 

IFIOABSlUJD.Gf  .OIFIGC  TO  49 

0061 

GO  TO  48 

0062 

49  Al-1. 

006  3" 

IF t  liJ  l.LT.O.  J  A1  =-l  • 

0064 

YIP=Y(  J-l  1*171  Jl-YI  J-ll  I*IU(l)*Al*OIF-U(J-l))/IUU)-UIJ-in 

0066 

GO  TO  44 

0066 

43  Y1P=0. 

c 

0067 

44  IFIXEX.NE.2)  GO  TO  45 

0068 

50  J=NP2 

0069 

51  J=J-1 

0070 

UJl-UI Jl-UINPiJ 

0071 

IF(0A8SIUJ1).GE.0(FIG0  TO  52 

0073 

52  Al-1. 

0074  ' 

IFIUJ l.LT.O. )A1=-1. 

0075 

YFM=Y<  J*lJ*tYlJ)-Y  I  J*l))*IUlNP3  J  +  Al*  DIF-UI.mil/IUIJI-UIJ4l>> 

00  76 

GO  TO  46 

0077 

45  YEM=YlNP3) 

0079 

RETURN 

0030 

END 

0001 

SUBROUTINE  M0MI2UMJ 

0002 

IMPLICIT  REAL*8I A-H»0— Z ) 

c 

M0HI2  CALCULATES'  INTEGRAL  MCMENTUM- TNCREM ENT~~IN~CHE6A  COORDINATES 

0003 

COMMON/V/  UI43) ,FI3 ,431,  RI43).  RH0I43),0MI43i, YI43J 

l/I/N,NPi;ht>2,Np3,KEQ,NPH,KEX,klM,<ASE,i<.kPAAN 


_ 0004  UD25=0.5*IUll)+UI3)-2 .*Ul NP31 ) 

0005 .  UDNP25=0.5*IUINPlI-tHNP3)  )  ' 

_ 0006  _ _  0  M2  5=0. 5*  OH  (3) 

0007  "  '  0MNP25=0.5*IQM(NP1  )♦!."»  ' 


C 

5TO5  Eht=0.5*0M25*tUD25*U( 1)-U(NP3I) 

0009  TO  10  I-4.NP1 

ooio  ioEHi-EHuo;r*roM(n-tinii-in*ium*oiT=rr=z^inNP3T) 

0011  EMI=EMI*0.5*0MNP25*U0AP25 

c  ‘  . . ' 

0012  AM=EMI 

0013  Return 


0014 


END 
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SUBROUTINE  OUTPUT'  " . 

IMPLICIT  REAL*8CA— HtO-Z J 

COMMON  /GEN/PE I, AMI  ,APE .DPOXi'PRE'FUT, PR  I 3)»P 1 3) . DEN, AMU. XU, XO.XP > 
lXL.OXfCSALEA.XPC6  >INTG 


)  tfl  3  IlKl.jltf  UItJll 

l/c/s: <43),AU(43> ,  BUI 43) ,CU(43 ), A{ 3,43  I .51 3,43 ) ,C I 3. 43) 
1/Ll /YL. UMAX. UMIN.FR.YIP.Y EM 

1/I/N,NP1.NP2.NP3.NEQ.NPH,KEX(KIN.KASE.KRA0.KPRAN 

C0MM0N/AUXP/TEMPI43) 

COMMON/ AUXY/YY 143 ) .XXU.RR1 


i  ;  ■va.i  t  1 1 1  :  1 1  v.  i  i  m  :  i  c*  1 1 


COMMON  /I DIN/  I NOIC 

C0MM0N/0U0/DU00MI43),  DUDY (43 ) .  A0UDY443).  AOUOYM 
COMMON  /ASD/  AS01.AS02 
"COMMON/ TEM/ TO jO 
COMMON/UMUM/YMU 


idiuniiMiLi*]  enimj  i 


COMMON/OUT /UP.US.PHIUCI 50 1) • PSIUI 43  I . X001 501 ) 

COMMON/TWOO/YI 

DIMENSION  AMACH 143 1 

GIMESrSlON  TEMPEI43) 

IF f INTG.NE.l)  GO  TO  15 


uu  11  1=1. 

AMACH II) =0. 

Tl  CONTINUE . 

WRITEI6.49) (0Mll),l=l,NP3> 

4?'F0RMATI24H1THE  VALUES  OF  CMEGA  ARE/I 11F10. 4)1 
15  CONTINUE 


UJU0*UI1)/UINP3) 

GO  TO  17 

16  UJUO*0. 

17  RHJO-RHOI 1I/RH0INP3) 
1FIKRAD.EQ.0I  GO  TO  605 


TEMPEII)-RII)/RO 
60  CONTINUE 
GO  TO  615 
605  DO  610  1*1, NP3 
610  TEMPEI I)*YII)/RO 


ntmiUiiaHii; 


AMACHIII-O. 

AMACH I NP3  >  *0. 

DO  62  1*2. NP2 

IFIOUDYI I I.EQ.O.I  GO  TO  62 

AMACHI  I)*SHEAR(II/IRHCII)*OUOYII)) 


HRI TE 16.51)  XXU.RR1 

51  FORMATI'O  XU*  ’.2P0U.2,'  RI  »  '.2P01L.2.’  IN') 

MRITF  16,551  UJUO ,RHJ  C.TOJO, PREF 1 1) , PREF (2)  .  PREFI3), AS01.AS02 
55  FORMAT  ( 1H0 , 

1  6HUJ/U0=F6. 3.2X.8HRHJ /RHC=F6 . 3 , 2X , 7HT0 J /T O* F6 .3 , 2X, 6HPREF1 


M  fVR MJ'IM  ^7 


wfMirMunmmtvii'vurax . 


23) 

HRITEI6.54) 

HRI TEI6.52) 

52  FORMA Tl  1H0.8X , 1HY ,11X,1HU.11X,1HH.12X, 1HC.9X.2HRE, 8X.4HR/R0 . 8X, 
15HSHEAR.9X.3HEPS.9X , 3HRH0 ,6X, 5HDU/0Y. 7X.4HPS IU/ 

29X.2HIN, 7X.6HFT/SEC.3X ,  12  HFT4*2/SEC**2» 

313X,11HIFT/SEC)**2 ,10X, 13HLBM/FT-SEC**2, 


EWffcH  l  A  A  £I^r^KWi:i  A  fi-1 afil 


53  FORMAT! 1H  1P11D12.3  ) 

54  FORMAT (1H0  ) 

00  10  J1-1.NP3 

J2=NP2-J1»2 

WRITE  1 6,53) YYIJ2) ,UIJ2),FI2,J2),FI3,J2),F( 1, J 2) , TEMPEl J2 ) , 


i  « 

10  CGNTINUE 
RETURN 
CND 


m  :  i  v  j  *  r*  rr  i :  li  * :  i  *  m  n  :  i  a  i  p  r4  n  ’  i  ’  i » i  i  p  m  it  it*  1 1 1 
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SUBROUTINE  RADIX, R1,CSALFBI 
IMPLICIT  REAL*B< A-H«C-2I 

APPLICABLE  TO  AXISYMMETFIC  MIXING  LAYER  AND  JET 

COMMON  /GEN/ PE  I ,AMI ,AME,OPDK, PREF (3 ) , PR  I  3) ,P 1 3) , PEN, AMU, XU, XO.XP, 
LXLiDXiCSALF A ,  X  P , 1 NTC 

1/V/UC43I ,F<3,4JJ,RI43I,FHC(43»,0MI43»,YC43I 
l/l/N,NPl,NP2,NP3,NE«,NPH,KEX,KIN,KASE,KRAD,  KPRAN 
COMMON/ RAOI/RO.RI 

common/thoo/yi 

commqn/iw/imake,isl 


IF  IKRAO.EQ.O)  GO  TO  18 
IFIKIN.E0.3)  GO  TO  17 
1FIX.EO.O.I  GO  TO  15 

Rl*R<l)*(ftm-2.*AMl*{X-XF>/IRHOm*U<ll  )  ) 
IFIR1.LT.0.IR1=0. 


18  Rl*l« 

IFIX.EQ.O.I  GO  TO  19 
IFIISL.EO.OI  RETURN 
IF( KIN.EQ.3)  GO  TO  20 
YI “YI-AHI *IX-XPI/IRH0I1I*UC1» I 


UI.LIa 

RETURN 

19  YI=RO-RI 
IF(ISL.EQ.U)  YI  *1. 
RFTURN 

20  YI=0. 
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SUBROUTINE  READY 
IMPLICIT  R£AL*8IA-H,0-ZI 

COMMON  /GEN/ PE  I .AMI , APE . OPDX , PREF (3)  ,  PR  (31  ,P  I  31 , OEN,  AMU,  XU.  XD  ,XP , 
1XL.0X.CSALFA.XPC6 ,INTG 

1/V/U(43),FI3,43I,RI43),RHCI43)«0M|43),Y(43)  ~ 

l/I/N.NPl, NP2 ,NP3,NE0,NPH,KEX,KIN,KAS E.KR AO, KPRAN 

common/thod/yi 

CALL  DfcNSTY 

CALL  RADIXU.RIil .CSALFAI 
C  Y  NFAR  THE  I  BOUNDARY 

~  '  IF  tPID.EQ.O.I  KlN=3  .  “ . .  ‘ 

IFIKRAO.Eg.il  GO  TO  70 
IFIYI.EO.O.)  KIN=3 

70  CONTINUE 
GO  TO  (71, 72, 731  , KIN 

71  CONTINUE 

GO  TO  74 - - - - -  -  '  - 

72  Yl2l=12.*0MI3)/( I3.*RHC(21*RH0(31 I* I U  12 1 +U I  3 1 +4.*UI 1) ) I 
GO  T  3  74 

73  YI2I*.5«CIM(3) /IRHO(l)  *U(11 1 

74  YI 31 *YI 21+. 25*OM( 31 *11./ I RHCI 31 *U 13)  1*2./ IRHOI 3»*Ut 31+RHOt 2)*UI  21  I 
II 

"C  "Y  *S"  FOR  INTERMEtTTATF  GRI0"-PCTNr5 -  - 

0018  00  50  1*4, NP1 

0019  50  YIII*YI(-1I*.5*I0MII1-CMI 1-1 11*  1 1 ./ IRHO  (I)  *U  I  1 1 1  *1.  / 1 RHO 1 1- 1 1  * 

1UI1— 1)11 

C  Y  NEAR  THE  E  BOUNDARY 

0020  Yl NP2 1 =Y I NPl I *.25*I OM INP2 1 -CM IN PI 1 1  *  1 1,/ I RHOI NP 1 1 *UI NP1 1  1*2./ 

•’  -  '  1IRH0INP1  l*UTNPll*RHOrNP2T*U(7iP2m -  - 

0021  GO  TO  1 81  ,82 • 83  I , KEX 

0022 
002  3 
0024 

”002? 

0026 
0027 
0029 
0029 
00  30 

'00  SI - 5TTW  54'  I^27nP3 . . .  ^ - 

0032  54  YIII=PFI*Y(I I /Rill 

0033  56  YI2)*2.*V(2)-YI3> 

C  CALCULATION  OF  RADII 
00  14  DO  57  I =2 »NP3 

0035  IFIKRAD.EQ.01RI I  1 =R 1 1  I 

00  36 - 1 F  ( KR  A  D  .  NE  .  CTTTTrn  ■TCm>Y  III  «CS  ALFA 

0037  57  CONTINUE 

0038  69  RETURN 

0039  END 


81  CONTINUE 
GO  TO  84 

82  Y|NP31=YINP2I*12.*ICMINP2 l-QHINPl )I/((RH0INP1|*3.*RHQ(NP2II*(U(NP2 
U*U(NPII*4.*UINP31I  I  . 

'  GC  TO  84“  -  -  - - 

83  YINP3l*YINP21*.5*I0M(NP2I-0M(NPll »/ IRHUINP3I*UINP3I I 

84  IF(CSALFA.£O.O..GR.KKAO.EQ.OI  GO  TO  51 
DO  52  1*2, NP3 

52  Y( 1 1 =2.*YC I) *PE1/ I R II l+USCRT I R 1 1 1 *R 1 1 1*2 .*Y 1 1 )*PE (*CSALFA 1 1 
GO  TO  56 


0001 

0002 

0003 


0004 
COO  5 
0006 

"0007  ' 

0008 
0009 
0010 
COIL 
0012 

0013 - 

0014 

0018 

0016 

0017 
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0001 

0002 

DOD'J 


0004 

DD05 

0DC6 

00U7 

0008 

COO? 

C010 

0011 

0012 

oou  — 

0014 

0015 

0016 

0017 

ooic 

~oonr  “ 

DD20 

0021 

0022 

0023 

0024 

“0025* . 

DO  26 
002  7 
0028 
0029 
0030 

-0031 - 

DO  32 
0033 
0034 
0035 
0D3& 

"OffST - 

0038 

C039 

0040 

0041 

0042 


0043 

D044 

D045 

0046 

0047 

0048 

0049 


SUBROUTINE  SHEARS 
IMPLICIT  RfAL*6(A-H,0-2) 

COMMON  /GEN/PEI .API .APE ,DPOX, PREF (31 , PR( 3 1 ,P 1 3). DEN, AMU, XU. XO , XP , 
1XL.0X. CSALFA .XPC6 ,  I  NTG 

1/ I/N, NP1 . NP2 ,NP3 , NEW, KPH.KEX ,KI N.KAS  b.KRAO, KPRAN 
1/V/UI43I »F(  3,43) ,R(43I ,  KHCI43I , CM (43  I ,V (431 
1/Ll /VL iUMAX.UMI N,FR,YIP,YEM 
C  CMMON/OUT/UP.US ,PHIUC(5CI).PS1U(43) . XOD ( 501 ) 

COMMON  /SHEAR/  SHF AR (431 .SCSHt 43 I 
COMMON  /ASD/  ASD1.ASD2 
C  OM MQN/AEMU/EMUA 
CCMM1N/0ENC /C 
C  CKMON/RT /RTM , YLK 

COMMON/OUO/OUODMI 43 1 ,  0UDYI43I,  AOUDY (43) ,  ADUDYH 
COMMON  /KUH /  RAAUHI43) 

C  OPMON/K JU/KMU 

.  C  CMMON/I  M/I  MAKE ,ISL 

COMMON/UMUM/YMU 
COMMON/ I 01 N/INDIC 
IKMAX=0 
TAUM-.OOOl 
AOUD YM= .0001 

-  00  98  1=2, NP1  . . 

RA=.5*(R( T*TI*R(I II 
RH=.5  *  I RHO I  I*1)+KHU(1 1 1 
UM>.5*(U(I»1)»U(I>) 

RAAUHl I ) =RA*RA*RH*UM 
SCSH(U*R(1)*RH0(II*U(II/PE1 

- DuoY(ir-ouoroMni*scsHnr-- - 

99  CONTINUE 

SCSH(  1)=R(  11*RH0(  1)*UI  1 l/PEI 

SCSHl NP2)=R (NP2 I *RHC( NP2)*U(NP2 l/PEI 
SCSHC  NP3I *R(NP3 )*RHOt  NP3) *U(NP3 l/PEI 
DUOYt  1 1 ■DUDOM (  1 1 *SCSH (  1) 

- DUmNP2,  =0UD0MTNP23  4SCSR1NP2I  - - 

0  UDV( NP3 1 =DUJQM( NP3 1  * SCSH ( NP3 I 
K  AAUH  ( NP2  I  =SCSH  (NP2I*MNP2)*PEI 
DO  97  1=1, NP3 
ADUOY (I  I =DA8S( DUOY ( 1 1 ) 

97  CONTINUE 

. . 00“  96  I*J4N . 

AOUDYL=AOUOYM 

IF  IADUOYI II.GT.ADUOYP)  AOUOYN- AOUDY (II 
IHADUDYH.GT.AOUDYLI  IMAX=I 
96  CONTINUE 

IF(INTfa.NE.l)  GO  TO  1066 

"  ~C‘ FI  171) ’INPUT  I  S’  EDDY'  VI5CCSTTT  ~FP5  -'STJ  TH  AT"TAU=RU(TFEFS*T30DT' 
f  ANO  TKE«EPS*0UDY/A1 
DO  1D45  I =2,NP2 

C  FOR  THIS  CALCULATION  INPUT  IS  EPS  (FT**2/SbC) 

1045  F(l,tl=F(lfl) *AOUDY (I l/ASDl 
F(1,1)=0. 

FTT7NP3)  so; - 

IF(K(N.EQ.3)F(1,1)=F( 1,21/2. 

XSTA=Y(NP3)*10. 

1CHK=7*NP3/10 
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0050 
0051 
00  52 
0053 
"  0  J54 
0055 
005& 

0057 

0058 

"  0009"~" 
OOoO 
004 1 
0062 
0063 
0064 
D065- 
0066 
0047 
0068 
0069 
0070 
0071" 
0072 
0073 
0074 
0075 
0076 

“6077 - 

0078 
0079 
0080 
0081 
0082 
""l)0"83’  _ 

00  P.4 
0085 
0086 
0087 
0088 
“0089 
0090 
0091 
0092 
0093 
0094 
"0095 
0096 
0097 
0098 
0099 
010C 
0101 
0102 
0103 
0104 


RTN=l. 

1066  CONTINUE 

IF! IHAX.GT. ICHKI  IMAX=ICHK 
SIGN=OUOVI 3 1 /OUDY ( NP1 )  _ 

DO  101  I “2«NP2 

33  SHEARm=4SDl*RH0II  »*F(1, 1  )*DUCY  II)/ AOUDV  II) 
(FIKlN.uQ.2l  GO  TO  101 
IF(  IWAKE.Ev).  II  GO  TR  333 
IF  I I.LT.IMAX.ANO.UIII .LT.0.9*UP) 

1SHC AK  ( I  )  =  SHE_AR  I 1  I  *ADUCYU  J/ADUDYM 
■  "GO  TO  lOf . . . .  . . . - 

333  IFISIGNIlOlti 34*334 

334  IFII.LT.i MAX]  SHE AR ( I )=SHEAR(I) *AOUOY ( 1 ) / ADUDYI IMAX) 

101  CONTINUE 
IFIKIN.EQ.3I  GO  Tn  118 
IFIX'J.LE.XSTAI  GO  TO  118 

- iPj urir."L'T.uTN"«)r  Gcr'rtrTD'z - 

IFIUHIN.LT. UINP3II  GO  TO  118 
GO  TO  103 

102  IFIUHIN.LT. Jill)  GG  TC  118 

103  CONTINUE 
DElU=U(NP3)-UU) 

rtf  i  I'd  f=i7W3  . 

IHU1  1I.GT.UINP3)  1  GO  TO  ill 
IFIIUII I-UI1I).GT..05*DELUI  GC  TO  112 
GO  TO  110 

111  IFI 1U1 I I-UINP31I.LT. -.95*DELUI  GO  TO  112 
110  CONTINUE 

112  t IN-I-I 

IF  I  1 1 N.LT .2 1  1 1 N=2 
DO  113  I«1,NP3 
IFIUI ll.GT.UI NP3I I  GO  TC  114 
IF|(UUI-U(1II.GT..95*DELUI  GO  TO  115 
GO  TO  113 

. in'  iFnum-uiNP3i».LT;--.054DEun  go  to  m 

113  CONTINUE 

115  IOUT=I 

IFI IPUT.EQ.NP3)  I0UT-KP2 

EPSl=SHEAR(I INI/I RHOI 1 1 Nl *DUDY( (INI  I 

F  PSO=SHE  AR  1 1  (JUT  )  /  (PHOI  ICUT)*DUDY  1 1  OUT  )  I 

00  116  "I “1«NM3 

IFI I.GT.I1N)  GO  TO  117 

SHEAR  ( I  )  =RH(J{  I )  *EPS  I  *DU0Y  1 1 1 

GO  TO  116 

117  IFI I.LT.IOUTI  GO  TO  116 
SHEAR ( I » =RHO( I ) *EP50*GUDY ( 1 1 

116  CONTINUE' 

118  CONTINUE 

Dtf  125  I “1«NP3 
TAUL=TAUM 

TAU=0A8S I  SHE ARI 1 1 1 
IFI TAU.GT.TAUH)  T  AUH=TAU 
I F  l  TA  UM.  GT.  T  AUL  )  I  KHAX  «I 
12S  CONTINUE 
RTHL=RTM 

C  HUH= SHEAR! I KHAX) / (RHCI IKMAX )*DUDY I IKHAX ) I 
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0105 
0106 
0107 
0108 
0109 
0110 
0111 
0112 
0113 
0114 
C  115 


0116 

0117 

one" - 

0119 

0120 

0121 

0122 

0123 

0124 - 

0125 

0126 

0127 

0128 


0129 

0130 

0131 

0132 

C133 


R  TM=DABS  I  III  1  >  -U  ( NP3J  l/EFUM 

YLK=1. 5708*DAHSCU(1)-U(NP3) )/ ADUDY ( I  KHAX ) 

IF( YLK.LT .YL)  YLK=YL 
1F(YLK.GT.1.57*YL)  YLK=1.57*YL 
1F1PSIU1I  KHAX) . LT.0.4 ) YLK=1. 57*YL 
1F1KIN.FU.3.AND.KRAD.NE.0)  YIK=2.*YMU 
P  TM=RTH*YLK 

IHXU.Lb.XSTA)  CO  TO  126 
IF! INTG.Eti.l)  KTML=RTH 
IF(RTM.GT.1.02*RTML>  RTM“1.02*RTML 
'  126 ‘CONTINUE  " 

C  ******* ************ *6** ****** ******** ********* 

C  RAMP  16  > 

C  *****************  ****** ****** ******** ********* 

IF(XU.LE.XSTA)  GO  TO  144 
I  FI R  TM.CT.185 • )  GO  TO  131 

-  A  $02=1.69  . . . 

GO  TO  145 

131  A  SD2* 0.464. 00 762 *RTP 

1F1AS02.GT.3.20)  ASD2=3.20 
GO  TO  145 

144  A  SD2=1 .69  ' 

-  145" CONTINUE  ‘ 

C  ********************************************** 

A  S02  =  ASD2/C 

S IGMA=l. 855*XU/YL 

WRITE (6 , 1000 11 NTG,RTM,ASD2, SIGMA 

1000  FORMAT  1 1H  FROM  SHEARSt  AT  STEP  ',14,'RTM*  '.1PE14.5 

- I1PE14.5;*'  SIGMA-  «  .IPE14.5T -  -  -  - 

WRI TF 1 6 ■ 1 001 )  YL.YLK.IKKAX 

1001  FORMAT!  1H  YL-  S1PE14.5,*  YLK-  S1PE14.5,'  UMAX-  *, 
FMUA=DABS (SHE AR 1 3 ) /DUDY 1 3 ) 1 
RETURN 
END 


•  ASD2-  '  > 
131 
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OOdl 

OOC2 

CJ03 


0004 

0005 


00U6 

0007 

0003 

0009 

0010 

0011 

0012 

0013 

0014 

0015 

0016 

0017 

CU3 

COl'i 

0020 

0021 

0022 

0023 

0024 

0025 

00?6 

0027 

0023 

0029 

0030 


0031 

C032 

0033 

0034 

0035 

0035 

0037 

0038 

0039 

0040 

0041 

0042 

0043 

0044 

0045 

0046 

0047 

0048 


SUBROUTINE  SLIP 
IMPLICIT  REAL*8(A-H,0-Z) 

(  (iMMMN  /GEN/ PEI  .AMI  ,  APE ,DPOX , PKEF (3 ) , PR (3 ) , P { 3) , DEN, AMU , XU,  XO,XP, 
.  1XL ,DX ,C  SALF  A ,XPC6 • I NTG 

1/ 1  /M,  NPl  •  NP2  ,NP3  ,NF  C,  NPH,  KfcX  ,KIN,  KAS  E , K.R  AO,  KP RAN 
1/V/llC  43)  ,F(3,43),R(43) ,RHC(43),0M(43),Y{43> 

COMMON/ACMU/ENUA 
COMMON  /L/AK.ALMG 

1/C/ SC (43) , AUI43) ,BU(43) ,CU( 43 ) , A I  3, 43 ) , B ( 3, 43 ) ,C ( 3, 43 ) 

C  SLIP  COEFFICIENTS  NEAR  THE  I  BOUNDARY  FOR  VELOCITY  EQUATION 
C  NOTE —  Ulll  MAY  OE  ZFRO  ONLY  IF  CPDX'l). 

CUI 2) =0. 

CU(NP2>*0. 

GO  TM  (71, 72,73), KIN 

71  CONTINUE 
GO  Til  74 

72  SQ=B4, *U( 1) *Ul 1)-12.*U( 1) *UI 3)»9.*U{ 3 )*U I  3) 
eu(2)=B. *(2.«U(ll+U(3))/(2.*U(l)»7.*U(i>*DSQRT(SQ>> 

AUI 2 I =1 .-BU( 2) 

GO  TO  74 

73  BUI  2 ) >0. 

1FIUI  D.EQ.O.)  GO  TC  735 
AKl=l./OX-0PUX/(RH0(l l*U( l)*U(l )) 

AK2=-U( 1 1 *AKl »DPDX/ (RHC (1 ) *U( 1) ) 

GO  TO  7355 
735  AK1=1./9X 
AK2*0. 

7355  CONTINUE 

AJ*RH1(1)*U(1)*,25*(V(2)»V|3) 1**2/ EMU  A 
IF(KPAU.CQ.O)  GO  TO  75 
AU(2)*2./(2.»AJ«AK1) 

CU(2)«-.5*AJ*AK24AU(2) 

GO  TO  74 

75  CU(2)~l./(2.»3.*AJ*AKl) 

AU(2)=CU(2) *(2.-AJ*AKl) 

C.U(2)=-CU(2)*4.*AJ*AK2 

C  SLIP  COEFFICIENTS  NEAP  THE  E  BOUND ARY  FOR  VELOCITY  EQUATION 
C  NOTF—  UINP3I  MAY  BE  ZERO  CNLY  IF  0PDX*0. 

74  GO  TO  ( 81 ,82  ,33) , KEX 

81  CONTINUF 
GO  TO  34 

82  S0*94.*U(NP3)*U(NP3I-12.*U(NP3)*U(NP1>»9.*U(NP1)*UINP1> 
AU(NP2)=8.*(2.*UI  NP3)  »U  (NPl )  ) /. { 2 .*U  (  NP3 )  »7.*U  (NP 1 J+DSQR  T(  SO  ) ) 
BUINP2I-1.-AUINP2) 

GO  TO  84 

83  AU(NP2)=0. 

CALL  V6FF|FMU,NP1,NP2) 

(F(U(NPi).CQ.O.)  GO  TO  835 

PK1*1« /UX-DPDX/ (RH0INP3) *U ( NP3) *U (NP3 ) ) 

BK2*~U( NP3) •3Kl»OPOX/ I RHC ( NP3 )*U (NP3 ) ) 

GO  TP  8355 
835  BKl=l./DX 
BKZ=*0. 

P355  CONTINUE 

EJ*RH0(NP3)*U(NP3)*.25*(2.*Y(NP3)-YINPU-Y(NP2))**2/EMU 

CU(NP2)-l./{2.»3.*BJ*BKl) 
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0049 

0060 

0051 

0052 
0053 
0054 
0055 
0056 
0057 
Ou63 
0059 
0060 
OCfal 
001  2 
0063 
00  64 
0065 
0066 
0067 
0068 
0067 
0070 
0071 
0072 
0073 
00  7  4 

0075 
0076 
0077 
C076 
0079 
oueo 
0091 
0082 
00  31 
0084 
0085 
0086 
0087 
0088 
0089 
0090 
0091 
0092 
0093 


8U(NP2)=CU(NP2)*t2.-QJ*BKll 
t U1NP2) *-CU 1NP2 ) *4. *8 J*BK2 
84  IF  1 NCQ.FQ. 1 1  RETURN 

C  SLIP  COEFFICIENTS  NFAR  THE  !  BOUNDARY  FCR  OTHER  EQUATIONS 
TO  54  J«1,NPH 
C  I J  .2 1  =0. 

C(J,NP21=0. 

CO  TO  (41,42.431 «KIN 

41  CONTINUE 
GO  TU  44 

42  A{J.2I-(U(2I*UI  3I-8.*UI  II  I /  I5.*(U<2  I *Ul  311+6.*U|  11) 

CF*l  l.-PREF  ( Jl  I  /  ( 1.4PRFM  Jl 1 
A(J,2l=<A(J,2]*r.FI/ll.4A(J,2l*GFI 

B(J,2!*l.-At J,2I 
GO  TO  44 

43  B ( J.2 I “0. 

call  sourceics.os.j.ii 

AKl-l./DX-DS 
AK2«-AK1*F( J.ll-CS 
AJF=AJ*PR£F (Jl 
(f IKRAD.EQ.Ol  GO  TO  45 
A( J,21=2./(2.*AJF *AK1  I 
C ( J t2 ) 5*A JF* AK2*A( J.2I 
GO  TO  44 

45  C( J.21*l./(2.»3.*AJF*AKll 
A(J,2)*CI  J(2I*I2.-AJF*AKI1 
C( J,21«-CI J,21*A.*AJF*AK2 

C  SLIP  COEFFICIENTS  NEAR  THE  E  60UNCAKY  FOR  OTHER  EQUATIONS 

44  GO  TO  151 (52 .531 .KEX 

51  CONTINUE 
GO  TO  54 

52  6 { J.NP21 * ( U( NP2 I ♦UINPll -8 .*U ( NP3 1 1/ lb.* IU INP2 I*UINP 1 1 1*8.*U(NP31 1 
GF=( 1 .-PREr I J 1 1 / 1 l.+PRE  F ( J) 1 

B(J,NP2I=(0(J,NP2 l+GF 1/ II. 401 J,NP2)*GFI 
AIJ.NP2l8l.-BIJ (NP21 
GO  TO  54 

53  AIJ.NP21-0. 

CALL  SJUKCb(CS(0S(J.NP31 

8K1-1. /OX-OS 

8K2«-BK1*F( J.NP31-CS 

ft  jF*fc  J*PREF  ( Jl 

C( J.NP2I ■l./(2.*3.*BJF*eKl| 

R  (  J  ,  NP2I  -C  (  J  ,  NP2I  *  (2.  — 13  JF  *BK1 1 
C I J . NP2 1 «-C IJ . NP2 I *4 . *B JF * BK2 

54  CONTINUE 
Ff TURN 
END 


0001 

0002 


0003 

0004 

0005 

0006 

0007 

TOf 

0009 

0010 

0011 

0012 

0013 

WIT 


SUBROUTINE  SOLVE (A.B.C.F, NP3I 
IMPLICIT  REAL *8 (A-H.O-Zl 
C  THIS  SOLVES  EQUATIONS  GF  THE  FORM 
C  Fill  -  A(I]*F(I*1)  *  B( 1 1 *F I 1— 1 1  «  C(l) 

X - FOR  I-2.NP2 - - - 

DIMENSION  AINP3I ,B (NP3I , C (NP31 , FINPJ 1 
NP2"NP3-1 

8(21  *  Bl  21  *F  (1 1  ♦  C 1 21 
00  48  I-3.NP2 
T  -  l./(l.-B(l) *A ( ] -1 1  I 

- A(  1 1  "*  "A  (I  J  *1 - 

48  8(11  -  IB((1*B(I-11  ♦  C 1 1 1  l*T 
00  50  1*2, NP2 
J*NP2-I»2 

50  FU1*AIJ1*F<J*1]>B{J1 
RETURN 
END 
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0001  SL6R0UTI NE  SOURCE  I CS ,DS >J , 1) 

0002  IMPLICIT  REAL*8I A-H*0-Z) 

C  FOR  CONSERVATION  OF  STAGNATICN  ENTHALPY 
C  CAUTION-  USE  CONSISTENT  UNITS 
~'T“  '  The 'DOT  PRODUCT  £}T  E  ViITH  4  IS  NEGLECTEO 

0003  CCMMON/GEN/PE I > AMI  ,  AME , GPOX , PREF I  3) . PR  I  3) ,P I  3 ) , OEN, AMU, XUtXDt XP , 

1XL.DX.CSALFA.XPC6, INTG 

1/V/UI 431 ,F<3,43),RI43),PHCC43),0MI43),YI43) 

1/I/N.NP1.NP2  ,NP3,NE<J,NPH,KEX,KlN,KASE,KRAD,KPRAN 


1/ LI /YL  >  UMA  X,UMIN»FK,Y[P,YEH 

T7C75CI43)  ,~AU  1431  7BU  Pf3 1 7C  UF4  3  J ,  A I  3,  4  3TrB73V4  31,  C I  3*431 


0004 

0005 

0006 

0007 

0008 

CCMMGN/AS0/ASU1 ,AS02 

COHMON/HT/RTN.YLN 

COMMON  /SHEAR/  SHEAR  1431 ,SC SHI 43) 

C0MM0N/DU0/0U0nMI43I ,  0UDYI43),  ADU0YI43).  AOUOYM 

COMMON/UMUM/YMU 

0009 

IF  ( JaGTai)  GO  TO  12 

0010 

GO  TO  (13,11, 12), J 

0011 

it 

IFII.EO.l)  GO  TO  12 

0012 

cs»scm*(ui  (*i)*u(i»i)-u(n*um>/toMiifri)-oMim 

0013 

cs=cs-sci  1-11  *iu<  n*um-uii-n*uu-i  n/ iomi  d-omi  i-ii  i 

0014 

CS-I 1. 0*0-1. O+O/PREFI  J) l+CS/ICMI l*l)-OMI 1-1)1 

“ 0TTT5 

rSKE=SCm'*  . '  TTTlVT*TI-RT7Tn/ 

1  IOMI 1*1 ) -UMI  I )  1 

0016 

CSKE*CSKE-SCI  I-1)*IF(  l,I)-FU,l-l))/(OMm-OM(  I-D) 

0017 

CS=CS+2.D*0*( l.U+O/PREF 111-1.0*0/ PREP IJ) )*CSKE/ IOM 1 I*1)-0MI 1-1) > 

0018 

CS=CS*CSKE 

0019 

OS-O.O+O 

Vl<  Mi  HP 

GO  TC)  3 

0021 

12 

CONTINUE 

0022 

CS  «  O.QD+O 

0023 

CS  *  O.OD+D 

"0024 

GO  TO  3 

C  SOURCE  MODIFICATION  NC.  1.  0«R*«2*RHO*U*EMU*DUDCM**2/PE I**2-DK/RHO*U 

13 

1FII.LT.2.DR. I.GT,NPZ7  GD  Tu  131 

0026 

CALL  VEFF (EMU, I ,1*1) 

0027 

GO  TO  132 

0028 

131 

CS=0.D*0 

0029 

GO  TO  133 

0030 

132 

CS-RHOIM  *U|I )*EMU  *I0U00MU)*R1I  )/PEI)**2 

133 

0032 

IFIUI 1) .EO. 0.0*0)  GO  TO  15 

0033 

C  S=C  S-ASD2+F  1 1 , 1)  **l.  50+0  / 1  YLK*U  ID) 

00  34  - 

DS*-1 . 5D+0*AS02*DSQRT (FC1 ,1)1/1 YLK*U II)) 

0035 

GO  TO  16 

0036 

- meat - 

15 

DS=0.D*0 

003  B 
0039 


16 


RFTURN 

ENO 


0001  SUBROUTINE  VEFF I  EMU , 1 , 1  PL  I 

0002  IMPLICIT  RE AL *8 1 A-H , 0-Z ) 

0003  COMMON  /SHEAR/  SHEAHI43) ,SCSH(43) 

0004  COMMON/OUO/OUOOMI43 1  *  0U0YI43),  ADUDYI43),  AOUDVM 

M05 - 67  Ehu4bABS(SHE*Rin/0UDYl  m - 

0006  RETURN 

0007  ENO 
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